Regression¶

Das Vorgehen in diesem Notebook orientiert sich am Data Science Lifecycle. Es wird eine Modellpräzision von > 70% gewünscht. Das bedeutet, dass der R-squared Wert des besten Modells >0,7 sein soll.

Python Setup¶

Zunächst werden mit der Import Methode wichtige Python Bibliotheken importiert.

import pandas as pd
import numpy as np 

import seaborn as sns
sns.set_theme(style='ticks')

import matplotlib.pyplot as plt
%matplotlib inline

import statsmodels.formula.api as smf
import statsmodels.api as sm
import statsmodels.gam.api as smg
from statsmodels.compat import lzip
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.tools.tools import add_constant
from patsy import dmatrix


from sklearn.compose import ColumnTransformer
from sklearn.compose import make_column_selector as selector
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.model_selection import train_test_split
import sklearn.metrics as skm
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Lasso
from sklearn.preprocessing import SplineTransformer
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline

import plotly.express as px
import plotly.graph_objects as go
import plotly.offline as pyo
# Set notebook mode to work in offline
pyo.init_notebook_mode()

Datenaufbereitung¶

Daten laden¶

Importieren der Daten aus der CSV in Github und Erstellung eines Pandas Dataframes zur Bearbeitung und Verwendung der Daten in diesem Notebook.

GIT = "https://raw.githubusercontent.com/jan-kirenz/project-annikas428/main/"
DATA = "project_data.csv"
TOKEN = "?token=GHSAT0AAAAAABPURMBOEW2WS3JEBXGI5TPAYPJRBIQ"

df_initial = pd.read_csv(GIT + DATA + TOKEN)

Daten inspizieren¶

Im ersten Schritt ist es wichtig, sich einen Überblick über die Daten zu verschaffen. Begonnen wird mit einer Anzeige einiger Datensätze.

df_initial

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	price_category
0	-122.23	37.88	41.0years	880	129.0	322	126	8.3252	452600.0$	NEAR BAY	above
1	-122.22	37.86	21.0	7099	1106.0	2401	1138	8.3014	358500.0	NEAR BAY	above
2	-122.24	37.85	52.0	1467	190.0	496	177	7.2574	352100.0	NEAR BAY	above
3	-122.25	37.85	52.0	1274	235.0	558	219	5.6431	341300.0	NEAR BAY	above
4	-122.25	37.85	52.0	1627	280.0	565	259	3.8462	342200.0	NEAR BAY	above
...	...	...	...	...	...	...	...	...	...	...	...
20635	-121.09	39.48	25.0	1665	374.0	845	330	1.5603	78100.0	INLAND	above
20636	-121.21	39.49	18.0	697	150.0	356	114	2.5568	77100.0	INLAND	above
20637	-121.22	39.43	17.0	2254	485.0	1007	433	1.7000	92300.0	INLAND	above
20638	-121.32	39.43	18.0	1860	409.0	741	349	1.8672	84700.0	INLAND	above
20639	-121.24	39.37	16.0	2785	616.0	1387	530	2.3886	89400.0	INLAND	above

20640 rows × 11 columns

Beim Sichten der Datensätze fällt auf, dass im ersten Datensatz zwei Unstimmigkeiten im Vergleich zu den anderen Datensätzen vorliegen:
median_house_value ist nur hier mit der Währung Dollar angegeben und hinter housing_median_age steht der Zusatz “years”.

Nun folgt Pandas Dataframe Info Methode. Diese zeigt die Variablen und deren Datentypen, sowie die Anzahl der non-null values je Spalte an.

df_initial.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 20640 entries, 0 to 20639
Data columns (total 11 columns):
 #   Column              Non-Null Count  Dtype  
---  ------              --------------  -----  
 0   longitude           20640 non-null  float64
 1   latitude            20640 non-null  float64
 2   housing_median_age  20640 non-null  object 
 3   total_rooms         20640 non-null  int64  
 4   total_bedrooms      20433 non-null  float64
 5   population          20640 non-null  int64  
 6   households          20640 non-null  int64  
 7   median_income       20640 non-null  float64
 8   median_house_value  20640 non-null  object 
 9   ocean_proximity     20640 non-null  object 
 10  price_category      20640 non-null  object 
dtypes: float64(4), int64(3), object(4)
memory usage: 1.7+ MB

Die Variablen housing_median_age und und median_house_value wurden als Datentyp object definiert. Dies muss genauer untersucht werden, da bei einer Jahreszahl bzw. einem monetären Wert ein numerischer Datentyp (int oder float) angemessener ist. Die Ursache hierfür ist auf die Zusatzangaben “$” und “years” zurückzuführen.
Die Variablen ocean_proximity und price_category sind keine numerischen sondern kategoriale Datentypen. Diese lassen sich jeweils in einer internen Rangordnung anordnen. Sie sind daher ordinal. Der richtige Datentyp hierfür ist category. Die Transformation der Datentypen erfolgt in Kapitel “Datentransformation”.
Neben den Datentypen wird auch die Anzahl der non-null values je Variable angezeigt. Hier haben alle Variablen eine Anzahl von 20640 non-null values. Lediglich die Variable total_bedrooms liegt unter diesem Wert. Dies lässt darauf schließen, dass hier einige Datensätze keinen Wert haben. Dies gilt es im Folgenden genauer zu untersuchen.

# show missing values -if present - in yellow
sns.heatmap(df_initial.isnull(), 
            yticklabels=False,
            cbar=False, 
            cmap='viridis');

print(df_initial.isnull().sum())

longitude               0
latitude                0
housing_median_age      0
total_rooms             0
total_bedrooms        207
population              0
households              0
median_income           0
median_house_value      0
ocean_proximity         0
price_category          0
dtype: int64

Wie bereits vermutet, enthält die Variable total_bedrooms 207 null values. Die Abbildung visualisiert diese nicht befüllten Datensätze und deren Verteilung im Datenset in gelb. Die null-values sind zufällig im Datenset verteilt. Da das Datenset über 20000 Datensätze enthält, werden die Zeilen mit den null-values im Kapitel Transform Data entfernt. Alternativ wäre auch eine Belegung der null-values beispielsweise mit dem Durchschnittswert der Variable möglich. Darauf wird hier jedoch verzichtet, genug andere Datensätze vorliegen.

Datentransformation¶

Entfernung der null values¶

# drop rows with missing values
df = df_initial.dropna()
print(df.isnull().sum())

longitude             0
latitude              0
housing_median_age    0
total_rooms           0
total_bedrooms        0
population            0
households            0
median_income         0
median_house_value    0
ocean_proximity       0
price_category        0
dtype: int64

Entfernung unwichtiger Features¶

Das Feature price_category wird entfernt, da es sich direkt aus der Variablen median_house_value ableitet und daher nicht zur Berechnung des Wertes genutzt werden darf. Dieses Feature wird im Notebook Classification als abhängige Variable genutzt.

#drop irrelevant features
df = df.drop(['price_category'], axis=1)

Konvertierung von Datentypen¶

Die mithilfe der info() Methode identifizierten Features, welche eine Datentyp Konvertierung benötigen, werden nun bearbeitet.

# Convert data types Ocean_Proximity and Price_Category to Data Type "Category" as they are categorical variables (ordinal)
df.ocean_proximity = df.ocean_proximity.astype('category')

# Housing_Median_Age and Median_House_Value are actually int64 values, but pandas did not get this because of signs like '€'
df.housing_median_age = df.housing_median_age.str.replace('years', '', regex=False)
df.median_house_value = df.median_house_value.str.replace('$', '', regex=False)
df.housing_median_age = pd.to_numeric(df.housing_median_age, errors='coerce').astype('int64')
df.median_house_value = pd.to_numeric(df.median_house_value, errors='coerce').astype('float64')
df.total_bedrooms = df.total_bedrooms.astype('int64')
df.info()
df.head()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 20433 entries, 0 to 20639
Data columns (total 10 columns):
 #   Column              Non-Null Count  Dtype   
---  ------              --------------  -----   
 0   longitude           20433 non-null  float64 
 1   latitude            20433 non-null  float64 
 2   housing_median_age  20433 non-null  int64   
 3   total_rooms         20433 non-null  int64   
 4   total_bedrooms      20433 non-null  int64   
 5   population          20433 non-null  int64   
 6   households          20433 non-null  int64   
 7   median_income       20433 non-null  float64 
 8   median_house_value  20433 non-null  float64 
 9   ocean_proximity     20433 non-null  category
dtypes: category(1), float64(4), int64(5)
memory usage: 1.6 MB

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity
0	-122.23	37.88	41	880	129	322	126	8.3252	452600.0	NEAR BAY
1	-122.22	37.86	21	7099	1106	2401	1138	8.3014	358500.0	NEAR BAY
2	-122.24	37.85	52	1467	190	496	177	7.2574	352100.0	NEAR BAY
3	-122.25	37.85	52	1274	235	558	219	5.6431	341300.0	NEAR BAY
4	-122.25	37.85	52	1627	280	565	259	3.8462	342200.0	NEAR BAY

Mithilfe der string replace Methode werden die zuvor identifizierten Zeichen “years” und “$” entfernt. Anschließend ist eine Transformation der Features in integer bzw. float Werte möglich. Für housing_median_age wird integer gewählt, da die Jahreszahlen als Ganze Zahlen also discrete Werte angegeben sind.
Für median_house_value wird der Datentyp float benötigt, da die Ergebnisse einer linearen Regression stehts continous sind.

Hinzufügen neuer Features¶

df = df.assign(people_per_household=lambda df: df.population / df.households)
df = df.assign(rooms_per_household=lambda df: df.total_rooms / df.households)
df = df.assign(bedrooms_per_people=lambda df: df.total_bedrooms / df.population)
df = df.assign(bedrooms_per_rooms=lambda df: df.total_bedrooms / df.total_rooms)

Diese Features werden dem Datenset hinzugefügt, in der Hoffnung, dass diese für die Verwendung in den Modellen nützlich sein könnten. Untersucht wird dies im Kapitel EDA.

Datensplit¶

Split des Datensets in 80% Trainingsdaten und 20% Testdaten. Die Testdaten werden anschließend bis zur finalen Evaluation der besten Modelle nicht verwendet.

train_dataset = df.sample(frac=0.8, random_state=0) 
test_dataset = df.drop(train_dataset.index)

train_dataset

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	people_per_household	rooms_per_household	bedrooms_per_people	bedrooms_per_rooms
14185	-117.08	32.70	37	2176	418	1301	375	2.8750	98900.0	NEAR OCEAN	3.469333	5.802667	0.321291	0.192096
6125	-117.91	34.11	20	3158	684	2396	713	3.5250	153000.0	<1H OCEAN	3.360449	4.429173	0.285476	0.216593
14095	-117.10	32.75	11	2393	726	1905	711	1.3448	91300.0	NEAR OCEAN	2.679325	3.365682	0.381102	0.303385
14359	-117.22	32.74	52	1260	202	555	209	7.2758	345200.0	NEAR OCEAN	2.655502	6.028708	0.363964	0.160317
18004	-121.99	37.29	32	2930	481	1336	481	6.4631	344100.0	<1H OCEAN	2.777547	6.091476	0.360030	0.164164
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
6168	-117.93	34.09	35	1891	353	1093	382	4.0167	165500.0	<1H OCEAN	2.861257	4.950262	0.322964	0.186674
4174	-118.21	34.10	36	2000	533	1234	535	3.7437	241700.0	<1H OCEAN	2.306542	3.738318	0.431929	0.266500
3593	-118.48	34.24	32	2621	412	1285	414	6.6537	267600.0	<1H OCEAN	3.103865	6.330918	0.320623	0.157192
5366	-118.39	34.04	44	1873	286	635	283	5.5951	461300.0	<1H OCEAN	2.243816	6.618375	0.450394	0.152696
17879	-122.00	37.40	17	5121	1017	1470	968	2.9706	81300.0	<1H OCEAN	1.518595	5.290289	0.691837	0.198594

16346 rows × 14 columns

Deskriptive Statistik¶

Kennzahlenanalyse numerischer Variablen¶

Mithilfe der describe() Methode werden die wichtigsten Kennzahlen numerischer Variablen (bspw. Durchschnitt, min und max Werte) ausgegeben.

#summary statistics for all numerical columns
train_dataset.describe().T

	count	mean	std	min	25%	50%	75%	max
longitude	16346.0	-119.564154	2.002618	-124.350000	-121.790000	-118.490000	-118.000000	-114.470000
latitude	16346.0	35.630318	2.138574	32.550000	33.930000	34.250000	37.710000	41.950000
housing_median_age	16346.0	28.664505	12.556764	1.000000	18.000000	29.000000	37.000000	52.000000
total_rooms	16346.0	2622.235776	2169.548287	11.000000	1448.000000	2119.000000	3120.750000	39320.000000
total_bedrooms	16346.0	535.281659	418.469078	3.000000	296.000000	432.500000	644.000000	6445.000000
population	16346.0	1416.087055	1103.842065	3.000000	784.250000	1164.000000	1711.000000	28566.000000
households	16346.0	496.758167	379.109535	3.000000	280.000000	408.000000	600.000000	6082.000000
median_income	16346.0	3.869337	1.902228	0.499900	2.555675	3.533200	4.744225	15.000100
median_house_value	16346.0	206916.154411	115676.394484	14999.000000	119300.000000	179700.000000	265900.000000	500001.000000
people_per_household	16346.0	3.089529	11.525259	0.692308	2.428850	2.816514	3.282584	1243.333333
rooms_per_household	16346.0	5.423747	2.160962	0.846154	4.443697	5.230149	6.051703	62.422222
bedrooms_per_people	16346.0	0.403370	0.223233	0.000670	0.315480	0.372646	0.443446	8.750000
bedrooms_per_rooms	16346.0	0.213188	0.057891	0.100000	0.175591	0.203263	0.239668	1.000000

Der Durchschnitt (mean) der Variablen housing_median_age über alle Datensätze beträgt 28,64 Jahre. Die Standardabweichung, also die durchschnittliche Abweichung der Werte vom Durchschnitt, beträgt 12,56 Jahre. Das Minimum liegt bei 1 Jahr, das Maximum bei 52 Jahren. 25% der Datensätze liegen unter 18 Jahren (1. Quartil). 50% unter 29 Jahren (2. Quartil) und 75% unter 37 Jahren (3. Quartil).
Im Durchschnitt beträgt der Wert people_per_household 3,09. Die durchschnittliche Abweichung von diesem Wert beträgt 11,53. Das Minimum liegt bei 0,69 und das Maximum bei 1243 Personen pro Haushalt. 25% der Datensätze liegen zwischen 0,69 und 2,43. 50% unter 2,82 und 75% der Datensätze unter 3,28 Personen pro Haushalt.

print(train_dataset.quantile(0.99, interpolation='nearest'))

longitude                 -116.290000
latitude                    40.630000
housing_median_age          52.000000
total_rooms              11021.000000
total_bedrooms            2208.000000
population                5723.000000
households                1955.000000
median_income               10.557500
median_house_value      500001.000000
people_per_household         5.286920
rooms_per_household         10.500000
bedrooms_per_people          0.941606
bedrooms_per_rooms           0.407472
Name: 0.99, dtype: float64

Die Quantile() Methode gibt für jede Variable einen Wert an, unter welchem (in diesem Beispiel) 99% der Datensätze liegen.
Beispielsweise liegen 99% der Datensätze unter einer Anzahl total_rooms in Höhe von 11021.
99% der Distrikte haben eine population unter 5723 Einwohnern.
Im Nachfolgenden werden einige mögliche Ausreißer identifiziert.

train_dataset.loc[train_dataset.people_per_household > 70]

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	people_per_household	rooms_per_household	bedrooms_per_people	bedrooms_per_rooms
9172	-118.59	34.47	5	538	98	8733	105	4.2391	154600.0	INLAND	83.171429	5.123810	0.011222	0.182156
19006	-121.98	38.32	45	19	5	7460	6	10.2264	137500.0	INLAND	1243.333333	3.166667	0.000670	0.263158
16669	-120.70	35.32	46	118	17	6532	13	4.2639	350000.0	NEAR OCEAN	502.461538	9.076923	0.002603	0.144068
3364	-120.51	40.41	36	36	8	4198	7	5.5179	67500.0	INLAND	599.714286	5.142857	0.001906	0.222222

Mehr als 70 Leute in einem Haushalten wirken sehr unrealistisch. 99% der Datensätze liegen unter einem Wert von 5,29 Personen pro Haushalt. Daher werden diese Datensätze als Ausreißer identifiziert.

train_dataset.loc[train_dataset.bedrooms_per_people > 8]

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	people_per_household	rooms_per_household	bedrooms_per_people	bedrooms_per_rooms
11862	-121.25	40.27	25	958	245	28	16	2.625	67500.0	INLAND	1.75	59.875	8.75	0.255741

Auch dieser Wert liegt weit über dem 99% Perzentil der Varialben bedrooms_per_people und wird daher als Ausreiser entfernt.

train_dataset.loc[train_dataset.rooms_per_household > 55]

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	people_per_household	rooms_per_household	bedrooms_per_people	bedrooms_per_rooms
12447	-114.49	33.97	17	2809	635	83	45	1.6154	87500.0	INLAND	1.844444	62.422222	7.650602	0.226059
11862	-121.25	40.27	25	958	245	28	16	2.6250	67500.0	INLAND	1.750000	59.875000	8.750000	0.255741
1913	-120.06	39.01	19	2967	528	112	48	4.0714	437500.0	INLAND	2.333333	61.812500	4.714286	0.177958

Auch diese Datensätze liegen weit über dem Wert des 99% Perzentils und werden daher aus dem Datenset entfernt.

#drop people per household outliers
train_dataset = train_dataset.drop([3364, 9172, 16669, 19006])
#drop bedroooms per people outliers
train_dataset = train_dataset.drop([11862])
#drop rooms per household_outliers
train_dataset = train_dataset.drop([1913, 12447])

#summary statistics for all numerical columns
train_dataset.describe().T

	count	mean	std	min	25%	50%	75%	max
longitude	16339.0	-119.564115	2.002470	-124.350000	-121.790000	-118.490000	-118.000000	-114.470000
latitude	16339.0	35.629562	2.138070	32.550000	33.930000	34.250000	37.710000	41.950000
housing_median_age	16339.0	28.664973	12.555983	1.000000	18.000000	29.000000	37.000000	52.000000
total_rooms	16339.0	2622.903544	2169.632035	11.000000	1449.000000	2120.000000	3121.500000	39320.000000
total_bedrooms	16339.0	535.416978	418.477283	3.000000	296.000000	433.000000	644.000000	6445.000000
population	16339.0	1415.032315	1100.485866	3.000000	785.000000	1164.000000	1711.000000	28566.000000
households	16339.0	496.956301	379.069221	3.000000	280.000000	408.000000	600.500000	6082.000000
median_income	16339.0	3.869002	1.901829	0.499900	2.555600	3.532900	4.744150	15.000100
median_house_value	16339.0	206925.109248	115665.635455	14999.000000	119300.000000	179700.000000	265900.000000	500001.000000
people_per_household	16339.0	2.941847	1.128881	0.692308	2.428863	2.816514	3.282400	63.750000
rooms_per_household	16339.0	5.413424	2.023751	0.846154	4.443662	5.230061	6.051152	52.848214
bedrooms_per_people	16339.0	0.402250	0.202972	0.014902	0.315552	0.372651	0.443430	7.925926
bedrooms_per_rooms	16339.0	0.213190	0.057897	0.100000	0.175587	0.203262	0.239667	1.000000

Analyse kategorialer Variablen¶

# summary statistics for all categorical columns
train_dataset.describe(include=['category']).transpose()

	count	unique	top	freq
ocean_proximity	16339	5	<1H OCEAN	7216

Die Variable ocean_proximity besitzt 5 Ausprägungen, wovon die Ausprägung <1H OCEAN mit einer Anzhl von 7216 am häufigsten vorkommt.

train_dataset['ocean_proximity'].value_counts()

<1H OCEAN     7216
INLAND        5215
NEAR OCEAN    2107
NEAR BAY      1796
ISLAND           5
Name: ocean_proximity, dtype: int64

Hier wird die Häufigkeit einzelner Ausprägungen von ocean_proximity dargestellt. Die Ausprägung ISLAND kommt über das gesamte Datenset mit über 20000 Datensätzen nur fünfmal vor. In Relation zur Gesamtanzahl der Datensätzen ist dies sehr wenig, wodurch die Aussagekraft dieser Ausprägung vermutlich eher gering ist. Daher werden die Datensätze mit Ausprägung ISLAND mit der Ausprägung NEAR OCEAN vereint.

train_dataset.loc[train_dataset.ocean_proximity == 'ISLAND']

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	people_per_household	rooms_per_household	bedrooms_per_people	bedrooms_per_rooms
8318	-118.48	33.43	29	716	214	422	173	2.6042	287500.0	ISLAND	2.439306	4.138728	0.507109	0.298883
8316	-118.32	33.33	52	2127	512	733	288	3.3906	300000.0	ISLAND	2.545139	7.385417	0.698499	0.240715
8315	-118.33	33.34	52	2359	591	1100	431	2.8333	414700.0	ISLAND	2.552204	5.473318	0.537273	0.250530
8314	-118.32	33.35	27	1675	521	744	331	2.1579	450000.0	ISLAND	2.247734	5.060423	0.700269	0.311045
8317	-118.32	33.34	52	996	264	341	160	2.7361	450000.0	ISLAND	2.131250	6.225000	0.774194	0.265060

# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR OCEAN': 'NEAR OCEAN',
    'NEAR BAY': 'NEAR BAY',
    'ISLAND': 'NEAR OCEAN',
    'INLAND': 'INLAND',
    '<1H OCEAN': '<1H OCEAN' }

train_dataset.ocean_proximity = train_dataset.ocean_proximity.map(remap_cat_dict).astype('category')

Exploratory Data Anaysis¶

Dieses Kapitel dient der ausführlichen Datenanalyse der Testdaten. Hierfür werden vor allem geeignete Visualisierungen in Form von Diagrammen verwendet.

erster Überblick¶

fig = go.Figure(data=go.Scattergeo(
        lon = train_dataset['longitude'],
        lat = train_dataset['latitude'],
        text = 'median house value in district: ' + train_dataset['median_house_value'].astype(str) + '$',
        marker_color = train_dataset['median_house_value'], marker_size=2, marker_colorscale = 'Inferno', marker_colorbar= {'title': 'Median House Value'}
        ))

fig.update_layout(
        title = 'Median House Value per District in California',
        geo_scope='usa',
    )

pyo.iplot(fig)

Die Datensätze stammen alle aus Kalifornien, USA.
Die Werte Longitude und Latitude unterscheiden sich daher kaum, wie sich auch bereits im Kapitel deskriptive Statistik gezeigt hat. Daher sind sie als Feature eher ungeeignet.
Allerdings scheint die Entfernung zum Ozean eine wichtige Rolle im Bezug auf den median_house_value zu spielen. ocean_proximity scheint daher als Feature interessant und sollte genauer untersucht werden.

train_dataset = train_dataset.drop(['longitude', 'latitude'], axis = 1)

sns.pairplot(train_dataset[['housing_median_age', 'total_rooms', 'total_bedrooms', 'population', 'households', 'median_house_value']])

<seaborn.axisgrid.PairGrid at 0x209e49bc850>

Error in callback <function flush_figures at 0x00000209E2DE5040> (for post_execute):

Die hier abgebildeten, unabhängigen Variablen weisen eine starke Korrelation untereinander auf (bis auf housing_median_age).
In Bezug auf die abhängige Variable median_house_value zeigen die Plots zunächst keinen starken Zusammenhang.
Interessant sieht das Histogramm von median_house_value aus, da es am äußeren rechten Rand nochmal hohe Häufigkeiten haben. Dieses Histogramm wird später nochmal genauer analysisiert.

sns.pairplot(train_dataset[['median_income', 'people_per_household', 'rooms_per_household', 'bedrooms_per_people', 'bedrooms_per_rooms', 'median_house_value']]);

Zwischen median_income und median_house_value ist ein linearer Zusammenhang erkennbar. Auch das Histogramm von median_income zeigt eine hohe Säule am rechten Rand.
Auch für die selbst definierten Variablen ist es schwierig, ein linearen Zusammenhang zu median_house_value zu erkennen.

fig = sns.displot(data= train_dataset, x='median_house_value')

Die Verteilung von Median_House_value ist nicht normalverteilt, rechts im hohen Preissegment ist die Häufigkeit nochmal extrem hoch. Über 500000 Dollar exisiteren keine weiteren Hauswerte. Möglicherweise wurde der Datensatz mit Werten über 500000 Dollar abgeschnitten und für alle höheren Werte, der Wert 500000 gesetzt.

fig = sns.displot(data= train_dataset, x='median_income')

Womöglich wurde auch die Variable median_income nach oben abgeschnitten.

Korrelationsanalyse¶

# Create correlation matrix for numerical variables
corr_matrix = train_dataset.corr()
corr_matrix['median_house_value'].sort_values(ascending=False)

median_house_value      1.000000
median_income           0.692099
rooms_per_household     0.176503
total_rooms             0.130828
housing_median_age      0.105534
bedrooms_per_people     0.077487
households              0.063096
total_bedrooms          0.047606
population             -0.026634
people_per_household   -0.178408
bedrooms_per_rooms     -0.253025
Name: median_house_value, dtype: float64

# Make a pretty heatmap

# Use a mask to plot only part of a matrix
mask = np.zeros_like(corr_matrix)
mask[np.triu_indices_from(mask)]= True

# Change size
plt.subplots(figsize=(11, 15))

# Build heatmap with additional options
heatmap = sns.heatmap(corr_matrix, 
                      mask = mask, 
                      square = True, 
                      linewidths = .5,
                      cmap = 'coolwarm',
                      cbar_kws = {'shrink': .6,
                                'ticks' : [-1, -.5, 0, 0.5, 1]},
                      vmin = -1,
                      vmax = 1,
                      annot = True,
                      annot_kws = {"size": 10})

Laut der Korrelationsmatrix einigen sich vor allem folgende Variablen:
median_income weist mit einem Wert von 0,69 die höchste Korrelation auf, danach folgen bedrooms_per_rooms mit einer negativen Korrelation in Höhe von -0,25 und anschließend people_per_household ebenfalls mit einer negativen Korrelation in Höhe von -0,18.
Danach folgt rooms_per_household, allerdings kann diese Variable nicht in Kombination mit bedrooms_per_rooms oder people_per_household genutzt werden, da sie sich aus total_rooms und households zusammensetzt. Diese Variablen werden bereits in bedrooms_per_rooms und people_per_household verwendet, daher ist die zusätzliche Verwendung von rooms_per_household nicht sinnvoll. Die Korrelation der Variablen untereinander ist in diesem Fall zu hoch mit einem Wert von -0,51 bzw. 0,6.
housing_median_age könnte mit einer Korrelation von 0,11 noch einen positiven Einfluss auf das Modell haben.

Analyse der kategorialen Variable Ocean_Proximity¶

sns.boxenplot(data=train_dataset, x='ocean_proximity', y='median_house_value')

<AxesSubplot:xlabel='ocean_proximity', ylabel='median_house_value'>

sns.stripplot(data=train_dataset, x='ocean_proximity', y='median_house_value', size=2)

<AxesSubplot:xlabel='ocean_proximity', ylabel='median_house_value'>

Die Ausprägungen NEAR BAY und NEAR OCEAN unterscheiden sich hinsichtlich median_house_value kaum. Daher können diese beiden Ausprägungen zusammengefasst werden.

sns.displot(data=train_dataset, x='median_house_value', hue='ocean_proximity', kind='kde')

<seaborn.axisgrid.FacetGrid at 0x230d8270f40>

Auch in der Häufigkeitsverteilung zeigen NEAR BAY und NEAR OCEAN kaum Unterschiede.

# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR OCEAN': 'NEAR COAST',
    'NEAR BAY': 'NEAR COAST',
    'INLAND': 'INLAND',
    '<1H OCEAN': '<1H OCEAN' }

train_dataset.ocean_proximity = train_dataset.ocean_proximity.map(remap_cat_dict).astype('category')

sns.boxplot(data=train_dataset, x='ocean_proximity', y='median_house_value')

<AxesSubplot:xlabel='ocean_proximity', ylabel='median_house_value'>

Auch NEAR COAST und <1H OCEAN zeigen eine ähnliche Verteilung des Medians und der Quartile. Das 3. Quartil von NEAR COAST zieht sich jedoch noch um ca. 50000$ weiter nach oben, als das von <1H OCEAN, daher bleibt die Differenzierung dieser beiden Ausprägungen noch bestehen.
Bei der Ausprägung INLAND sind viele Ausreiser nach oben zu erkennen.

Analyse numerischer Variablen¶

# check relationship with a joint plot
sns.jointplot(x="median_income", y="median_house_value", hue="ocean_proximity", data=train_dataset);

# check relationship with a joint plot
sns.lmplot(x="median_income", y="median_house_value", data=train_dataset);

Wie bereits aus der Korrelationsmatrix zu schließen, ist ein linearer Zusammenhang zwischen median_income und median_house_value erkennbar. Die “Linie” mit Datenpunkten bei 500000 Dollar ergibt sich aus der anzunehmenden Begrenzung der median_house_value_ auf 500000 Dollar.
Auch innerhalb der einzelenden Ausprägungen von ocean_proximity ist eine lineare Korrelation von median_income und median_house_value erkennbar.

sns.relplot(data=train_dataset, x='total_rooms', y='median_house_value', hue='ocean_proximity', col='ocean_proximity')

<seaborn.axisgrid.FacetGrid at 0x193868f47c0>

Ein linearer Zusammenhang von total_rooms zu median_house_value ist nicht wirklich zu erkennen. Die Variable ist daher für eine lineare Regression eher ungeeignet.

# Plot with Plotly Express
fig = px.scatter(x=train_dataset['bedrooms_per_rooms'], y=train_dataset['median_house_value'], opacity=0.65, 
                trendline='ols', trendline_color_override='darkred')
pyo.iplot(fig)

fig = px.scatter(x=train_dataset['people_per_household'], y=train_dataset['median_house_value'], opacity=0.65, 
                trendline='ols', trendline_color_override='darkred')
pyo.iplot(fig)

bedrooms_per_rooms und people_per_household zeigen einen negativen linearen Zusammenhang. Bereits in der Korrelationsmatrix war der negative Zusammenhang dieser beiden Variablen zu median_house_value sichtbar. Die Variablen werden nachfolgend in die linearen Regressionsmodelle aufgenommen.

Feature Auswahl¶

Auf Basis der in der deskriptiven Statistik und EDA gewonnen Erkenntnisse, sollen nun geeignete Features für die Modellierung ausgewählt werden.
Diese Auswahl wird auf Basis der sogenannten Filter Methode getroffen.
Diese eine von drei Möglichkeiten bei der Feature Selektion. Dabei werden beispielsweise auf Basis einer Korrelationsmatrix die besten Feature ausgewählt und nur diese dem Modell zur Verfügung gestelltn.
Folgende Features wurden als geeignet identifiziert:

median_income
ocean_proximity
bedrooms_per_rooms
people_per_household
housing_median_age (eventuell)

Modellierung¶

Statsmodel¶

Model 1 - OLS Regression¶

Als erstes Modell wird statsmodels OLS Regression verwendet.
Im ersten Modell werden lediglich die Feature 1-4 verwendet.

# Fit Model
lm = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household' , data=train_dataset).fit()
# Full summary
lm.summary()

OLS Regression Results
Dep. Variable:	median_house_value	R-squared:	0.623
Model:	OLS	Adj. R-squared:	0.623
Method:	Least Squares	F-statistic:	5407.
Date:	Wed, 19 Jan 2022	Prob (F-statistic):	0.00
Time:	12:34:59	Log-Likelihood:	-2.0569e+05
No. Observations:	16339	AIC:	4.114e+05
Df Residuals:	16333	BIC:	4.114e+05
Df Model:	5
Covariance Type:	nonrobust

	coef	std err	t	P>\|t\|	[0.025	0.975]
Intercept	3.394e+04	4569.346	7.428	0.000	2.5e+04	4.29e+04
ocean_proximity[T.INLAND]	-6.785e+04	1408.255	-48.179	0.000	-7.06e+04	-6.51e+04
ocean_proximity[T.NEAR COAST]	1.557e+04	1423.620	10.934	0.000	1.28e+04	1.84e+04
median_income	4.323e+04	404.027	107.007	0.000	4.24e+04	4.4e+04
bedrooms_per_rooms	3.009e+05	1.3e+04	23.197	0.000	2.75e+05	3.26e+05
people_per_household	-1.376e+04	496.667	-27.713	0.000	-1.47e+04	-1.28e+04

Omnibus:	5010.041	Durbin-Watson:	1.972
Prob(Omnibus):	0.000	Jarque-Bera (JB):	27166.062
Skew:	1.373	Prob(JB):	0.00
Kurtosis:	8.689	Cond. No.	128.

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

lm.params

Intercept                         33939.853081
ocean_proximity[T.INLAND]        -67848.787400
ocean_proximity[T.NEAR COAST]     15566.550456
median_income                     43233.585628
bedrooms_per_rooms               300850.020547
people_per_household             -13763.891665
dtype: float64

Die summary() Methode berechnet die wichtigsten Kennzahlen des Modells.
R-squared ist eine Kennzahl zur Beurteilung der Anpassungsgüte eines Modells und nimmt einen Wert zwischen 0% und 100% an. Als Bezugsbasis wird der Durchschnitt verwendet.
lm erreicht einen R-squared Wert in Höhe von 62,3%. Somit können 62,3% der median_house_value Streuung
In statistischen Methoden wird meist lieber der adjusted R-squared genutzt. Da dieser die degrees of freedom miteinbezieht, lässt sich durch adj. R-squared einen besseren Rückschluss auf die Gesamtpopulation der Datensätze ziehen. Dieser ist meistens etwas schlechter als R-squared. In diesem Modell sind beide Werte gleich hoch.
Mithilfe der Kennzahlen AIC und BIC soll ein overfitting des Modells verhindert werden. Für jedes zusätzliche Feature wird eine “Bestrafung” in die Kennzahlen miteinberechnet. Je kleiner die Werte, desto besser. In diesem Modell sind die Werte jedoch relativ hoch.
F-statistics ist die Menge der systematischen Varianz (MSM) geteilt durch die Menge der unsystematischen Varianz (MSR). Die Kennzahl gibt an, in welcher Höhe das Modell die Ergebnisausgabe der abhängigen Variable verbessert hat, verglichen zu Ungenauigkeit im Modell. Je höher F-statistics, desto besser. Ein F-statistics Wert in Höhe von 5407, wie in Modell lm, ist schon sehr gut.
Der Y-Achsenschnittpunkt des Modells liegt bei 33939,85. Die Steigung der einzelnen Features ist in der Spalte coef abzulesen. Zu interpretieren sind diese folgendermaßen: Vergrößert sich der Wert in Variable median_income um 1 und alle anderen Werte bleiben konstant, so ist das Ergebnis um 43233,59 größer als zuvor.
Die p-values der Features liegen alle bei 0,00. Die Nullhypothese (ß=0 => Feature hat keinen Einfluss auf die abhängige Variable) kann demnach für alle Features verworfen werden.
Die p-values für Jarque-Bera und Omnibus liegen ebenfalls bei 0,00. Die Nullhypothese, die Fehler seien normalverteilt, wird in diesem Falle verworfen. Es wird also angenommen, die Residuals seien nicht normalverteilt. Bei einer Datensetgröße > 50, kann dies aufgrund des Central Limit Theorems jedoch vernachlässigt werden.
Der Durbin-Watson Wert liegt bei 1,97. Somit kann angenommen werden, dass die Fehler des Modells keine Korrelation aufweisen.

Regression Diagnostics¶

fig = sm.graphics.influence_plot(lm, criterion="cooks")
fig.tight_layout(pad=1.0)

Mithilfe des Influence Plots werden mögliche Ausreißer als blaue Punkte in der Größe ihres Einflusses identifiziert. Dies können Werte mit hoher Abweichung in ihrem Y-Wert (outliers) oder X-Wert (high leverage) sein. In einer linearen Regression können solche Ausreise einen großen Einfluss auf das Gesamtmodell haben und somit die Performance des gesamten Modells erheblich verschlechtern. Als Faustregel werden Datenpunkte mit einer Cooks Distance in Höhe von 4/n (wobei n die Anzahl der Datensätze ist) als mögliche Ausreißer identifiziert.

# obtain Cook's distance 
lm_cooksd = lm.get_influence().cooks_distance[0]

# get length of df to obtain n
n = len(train_dataset["median_income"])
print('n = ', n)

# calculate critical d
critical_d = 4/n
print('Critical Cooks distance:', critical_d)

# identification of potential outliers with leverage
out_d = lm_cooksd > critical_d

# output potential outliers with leverage
print(train_dataset.index[out_d], "\n", 
    lm_cooksd[out_d])

n =  16339
Critical Cooks distance: 0.0002448130240528796
Int64Index([15778,  6776, 11511,  7263, 15723, 14756, 15687,  4638, 14806,
            15815,
            ...
            15809,  8866,   199, 18498,    59, 10749,  8857, 13154,  5470,
            16642],
           dtype='int64', length=934) 
 [1.49490398e-03 2.68414884e-04 4.78133170e-04 5.02665948e-04
94624792e-04 1.45661981e-02 6.06459084e-04 9.89584645e-04
26161610e-04 8.30878122e-04 1.36260242e-03 9.64095843e-04
86783348e-04 3.28695576e-04 3.75058682e-04 7.93946113e-04
82170067e-04 2.59400713e-04 5.38430141e-04 4.34597064e-04
58773707e-04 4.39852205e-04 3.24124357e-04 3.81088866e-04
95276469e-04 8.44927027e-04 4.83475406e-02 5.21391591e-04
36852748e-04 4.06798688e-04 3.74228212e-04 2.75462275e-04
37885379e-04 4.91036338e-04 1.65550557e-03 2.28062263e-03
11362268e-03 5.40836931e-04 5.06680630e-04 2.80028526e-03
84825684e-04 4.88676566e-04 8.68935319e-04 2.51788726e-03
76056994e-04 3.19987726e-04 3.94500298e-04 8.87693727e-04
88788562e-04 4.33792026e-04 3.43876382e-04 4.28843898e-04
54230456e-04 3.50021703e-04 4.80369640e-04 5.30319765e-04
59051144e-04 6.00454845e-04 3.80966666e-04 4.31199706e-04
03345609e-04 8.57174195e-04 4.11096007e-04 4.63221951e-04
23566980e-03 2.50497036e-04 4.82794226e-04 5.26525139e-04
69167729e-04 3.28935485e-04 6.97236942e-04 1.22926028e-03
37881318e-04 2.93996391e-04 3.92057182e-04 4.70614674e-04
17905418e-03 5.23060291e-04 4.86889004e-04 5.40614399e-04
95440918e-03 2.90055671e-04 9.44897532e-04 7.33532196e-04
47125775e-04 3.18823752e-04 1.06578506e-03 3.58388250e-04
73670560e-04 1.10228431e-03 7.93461030e-04 2.61856821e-04
98663391e-04 5.60713494e-04 8.96529290e-04 2.49653913e-04
07462907e-03 3.66163173e-04 3.27982973e-04 6.74318828e-04
07600427e-03 3.12835081e-04 9.41739908e-04 1.06035239e-03
24403292e-04 2.32717050e-03 3.33725982e-04 3.31713778e-04
45370019e-04 5.19796814e-04 2.52998567e-04 3.15019545e-04
17730666e-04 4.95060839e-04 5.84604429e-04 6.76035191e-04
62979595e-04 2.78405274e-04 2.64614306e-04 3.43991908e-04
28547805e-04 3.96502354e-04 5.72061896e-04 9.03304907e-04
88962274e-04 3.88578207e-04 6.03733226e-04 4.42201426e-04
84262789e-04 3.02906624e-04 6.67400127e-04 4.21938365e-04
18350166e-02 2.50731877e-04 7.56862933e-04 3.26525900e-04
91886314e-04 4.94062027e-04 2.62415163e-04 7.40021481e-04
74038819e-04 6.64113946e-04 6.00653167e-04 2.59384025e-04
68021977e-04 5.08871869e-04 3.45793984e-04 3.16681659e-04
59376409e-04 2.75726762e-04 5.71869366e-04 2.55386427e-04
78580587e-04 4.80446565e-04 6.29032892e-04 2.68844084e-04
15777697e-04 4.18502209e-04 2.66244188e-04 9.72580016e-04
59936346e-04 1.26716166e-03 4.61208595e-03 2.72426608e-04
31337333e-04 6.13153738e-04 3.94343301e-04 3.35424565e-04
11719688e-04 2.87763886e-04 2.48788511e-04 5.29875909e-04
52462604e-04 6.18585440e-04 2.93302623e-04 2.59391598e-04
57900725e-04 2.67878340e-04 7.54047576e-04 7.90100598e-04
51351563e-04 8.86252032e-04 1.18258562e-03 2.64701347e-04
06120700e-03 5.21490540e-04 3.94491160e-04 2.69647319e-04
65030920e-04 1.24489084e-03 6.75547522e-04 3.70755646e-03
64789182e-04 2.46399968e-04 3.96264268e-04 2.99146468e-04
65433172e-04 2.98636628e-03 3.31239647e-04 3.41387971e-04
69398818e-04 6.40829474e-04 3.15785273e-03 2.61734223e-04
60648977e-04 2.75598565e-04 4.26681256e-04 1.74170377e-02
74987588e-03 2.79744243e-03 8.47589712e-04 1.34331315e-03
85411722e-04 4.38351271e-04 3.92230730e-04 2.69865195e-04
13534479e-04 6.00850811e-04 2.56881788e-04 3.00847558e-04
60229616e-04 3.71622289e-04 4.01045487e-04 4.40508578e-04
71587095e-04 3.40572615e-04 2.68812469e-04 3.70078254e-04
10179623e-04 3.02152403e-04 2.67996019e-04 7.44260318e-04
13479879e-04 8.19067610e-02 5.12803460e-04 3.93272586e-04
19823559e-03 3.53653351e-04 4.34626944e-04 1.39726947e-03
18617182e-03 5.01736586e-04 5.30094203e-04 3.33634176e-04
18932862e-02 9.42177142e-04 4.18773467e-04 6.03680488e-04
71553393e-04 3.74286350e-04 1.23669485e-03 4.65135400e-04
69268542e-04 4.13374391e-04 3.53588165e-04 9.81277198e-04
54137328e-04 1.35995340e-03 3.11246270e-03 4.31574146e-04
86358637e-04 2.69481369e-04 2.98644303e-04 5.35187123e-04
51539827e-04 3.13073701e-04 9.94563742e-04 3.79299795e-04
72704600e-04 2.61286323e-04 2.94840250e-04 2.61668973e-03
65733675e-04 3.11871682e-04 1.16617792e-03 3.22950323e-03
59684839e-04 1.61885963e-03 4.58416247e-04 2.48030609e-04
42383137e-04 4.48068631e-04 6.04910702e-04 8.04956576e-04
11900750e-04 2.50644499e-04 3.76788326e-04 5.23889548e-04
98319955e-04 2.47144472e-04 3.04257523e-03 2.97617992e-04
70824454e-04 5.76938158e-04 7.24114626e-04 3.09400255e-03
36686516e-03 4.42145065e-04 3.02082029e-04 2.63000798e-04
88936728e-04 3.38298305e-04 4.44264842e-04 2.71776042e-04
56975665e-04 4.40378588e-04 9.46502316e-04 4.44336623e-04
68241357e-03 8.19569279e-04 7.25573623e-03 4.63097728e-04
86542144e-03 1.79765644e-03 3.13796286e-04 1.85849209e-03
63517011e-04 3.05689630e-04 4.46005083e-04 1.86968485e-03
21898380e-04 6.28339395e-04 2.74402091e-04 5.88073246e-04
32709199e-03 4.52837607e-04 1.50790793e-03 5.81240516e-04
64359329e-03 5.80674834e-03 6.55261768e-04 1.73632609e-03
32072802e-04 4.10811940e-04 1.04001291e-03 5.98785387e-04
79515437e-04 6.23027502e-04 4.23275012e-04 1.38618438e-03
71617878e-04 2.39674723e-03 3.51234700e-04 3.86146581e-04
07595499e-04 4.20280178e-03 3.69799077e-04 2.58852189e-04
00625914e-03 2.14059605e-03 5.24085589e-04 5.86927033e-04
76194574e-04 6.77176762e-04 6.22637307e-04 6.22956497e-04
12480619e-03 4.16146778e-04 3.70163512e-04 7.00008802e-04
59599089e-04 3.82851777e-04 4.55744907e-04 1.10934093e-03
52817097e-04 2.87570210e-04 2.52329760e-04 7.02649092e-04
59392154e-04 4.46540428e-04 5.57525212e-04 4.63626144e-04
65869117e-04 5.78289495e-03 6.15022887e-04 3.35508687e-04
65061093e-04 2.02798819e-03 2.49284001e-04 3.23955147e-04
10369090e-04 3.80530786e-04 4.15664656e-04 2.95094039e-04
47132464e-04 2.85414761e-04 2.84178883e-04 2.65532294e-04
71743672e-04 3.20261144e-04 4.67532383e-04 3.33304384e-04
52355395e-04 1.28385469e-03 2.87227641e-04 7.19836520e-04
24407306e-03 3.60495943e-04 2.99981365e-04 4.65103270e-04
46220633e-04 4.25385196e-04 1.69057212e-03 9.95180258e-04
27190396e-04 2.85787247e-04 3.28116249e-04 7.50507763e-04
31914472e-04 3.33314939e-04 3.96174751e-04 4.33313938e-04
35842338e-04 3.91664204e-04 5.53686911e-04 4.22118493e-04
49932882e-03 7.81329389e-04 2.75321278e-03 2.71498687e-04
92285936e-04 1.11281331e-03 2.80607976e-04 4.96095023e-03
21027426e-04 3.23368164e-04 4.14159128e-04 3.77467629e-03
07651838e-04 7.18666140e-04 4.53756366e-04 3.41057171e-03
77265109e-04 2.55933137e-03 3.99882274e-04 5.34894447e-04
20700768e-04 2.52589779e-04 8.76764590e-04 6.01984749e-04
01172240e-03 4.34738639e-04 4.86857764e-03 4.44271558e-04
77835060e-04 4.17273122e-04 2.49841045e-04 8.57922366e-04
01138406e-03 3.17770892e-04 2.79425670e-04 2.87924766e-04
75321438e-04 2.16673451e-03 5.56504395e-04 3.85884793e-03
33212441e-04 6.84666057e-04 3.21755993e-04 5.23263229e-04
07589874e-04 1.36587639e-03 8.41514470e-04 6.10850040e-04
87242598e-04 3.28588051e-01 3.71667189e-04 2.79196718e-03
09579459e-04 4.37807917e-04 8.51212481e-04 4.86433645e-04
30857717e-04 5.01044984e-04 2.47930811e-04 2.93683726e-03
59552309e-04 2.57442888e-04 1.45592560e-03 3.45294621e-04
20625186e-04 2.50894360e-04 7.58404371e-04 3.63506552e-04
32150175e-03 3.57160051e-04 2.82166230e-04 4.05308960e-04
05660916e-03 2.36117472e-03 2.55300114e-04 2.99426708e-04
48755421e-04 3.47614085e-04 5.35389815e-04 2.28333284e-03
69689693e-04 1.83993869e-02 2.78134600e-04 2.94646617e-04
32580902e-04 3.37268213e-04 8.49437828e-04 4.70592108e-04
65122577e-04 7.87353472e-04 2.65151373e-04 3.40147104e-04
47373921e-03 3.58812729e-04 6.71821996e-03 3.17543351e-04
61536099e-04 9.04605585e-03 5.89081995e-04 7.59534268e-04
64907542e-04 2.85244797e-03 6.43833067e-03 3.42696744e-04
53282313e-04 3.02325255e-04 8.96809026e-04 2.57985742e-04
29591390e-04 2.56536173e-03 3.98997462e-04 3.80381504e-04
04846515e-04 5.53919005e-04 1.72091124e-03 3.41707784e-04
54876276e-04 8.50177014e-04 3.12956299e-04 1.32323097e-03
46349028e-04 3.41588528e-04 5.65716565e-04 5.51966438e-04
49538242e-04 5.01573516e-04 2.67487102e-04 7.86724676e-04
66813565e-04 2.81069561e-04 3.72920665e-04 2.67001014e-04
45976115e-04 8.53312926e-04 3.47462338e-04 3.12518224e-04
37634912e-04 6.71193185e-03 2.80483340e-04 1.05535803e-03
05727691e-04 5.49198981e-04 3.92816712e-04 3.26462785e-03
78658807e-04 1.07133515e-03 2.58061845e-04 6.28308148e-04
99432217e-04 4.08464699e-04 7.76299893e-03 3.32998562e-04
75864814e-04 2.43445894e-03 1.29339473e-03 3.88447003e-04
72155160e-04 1.16799456e-03 6.99554601e-04 3.41654550e-04
60497617e-04 6.14377217e-04 1.12340141e-03 9.52380819e-04
28633882e-04 2.66993087e-03 7.63434605e-04 2.72831557e-04
40328575e-04 3.27503237e-04 9.76400663e-04 8.92877629e-04
15770615e-04 3.34734715e-04 6.30237071e-04 3.76652904e-04
53213140e-04 4.80060218e-04 3.10904734e-04 5.28275910e-04
94585293e-04 2.14279533e-03 3.85897937e-04 5.67212170e-04
67990920e-04 4.60165209e-04 7.28971845e-04 8.93621528e-04
40300668e-04 5.67891373e-04 2.98114027e-03 5.06945094e-04
90940363e-04 7.08720513e-04 2.58135066e-04 2.93029122e-04
81941867e-04 9.58097289e-04 5.82149611e-04 4.05186998e-04
70767144e-03 2.78477079e-03 2.89663179e-04 1.20228109e-03
54930522e-04 2.70933886e-04 4.04780848e-04 7.11395439e-04
32404421e-04 2.79054479e-04 3.51763013e-03 3.66814342e-04
17607727e-04 1.92440482e-03 7.40987676e-04 6.02952970e-04
90803016e-04 2.87937670e-04 4.13298793e-04 6.02532557e-04
03667841e-04 2.85359319e-04 4.43822131e-04 4.48611300e-04
84253469e-03 2.45013962e-04 5.37653017e-04 3.52961018e-04
78910376e-04 1.02483048e-03 1.04860911e-03 1.49712177e-03
69398997e-04 4.40503340e-04 2.55594132e-04 1.19645327e-03
53989647e-04 4.47037845e-04 3.51513607e-04 4.08201711e-04
50074578e-03 8.30344714e+00 5.65625002e-04 2.46168492e-04
33248865e-04 1.22758173e-03 1.01208124e-03 4.25612245e-04
64096985e-04 2.47406741e-04 2.70997639e-04 2.97073012e-04
13625041e-04 1.40881899e-03 3.02927343e-04 7.98328272e-04
70190253e-04 3.24999226e-04 3.06698993e-04 3.93282175e-04
77023978e-04 6.30043037e-04 9.15404010e-04 3.07118690e-04
75683525e-04 2.69553212e-04 6.54946966e-04 2.95261234e-04
85318203e-03 5.86215552e-04 3.24299295e-04 9.44380999e-04
08059729e-04 1.78452639e+00 2.78842121e-03 2.93019660e-04
81763831e-04 2.46515742e-04 2.67421853e-04 2.46015947e-04
59353083e-04 3.30689996e-04 3.10118580e-04 5.66627749e-04
76468283e-04 5.23602207e-04 3.17035827e-04 4.32554132e-04
88423909e-03 3.18177800e-04 2.47097467e-04 8.97296353e-04
11799006e-03 3.76872367e-04 5.61415054e-04 4.44857116e-04
99348646e-03 2.51524523e-04 5.19464921e-04 5.07202394e-04
18908548e-03 4.46062016e-04 1.23169781e-03 4.41564831e-04
31238292e-04 3.78250481e-04 2.59072590e-04 2.57755138e-04
56018898e-04 5.76159655e-04 1.67865668e-03 2.98192987e-04
79505548e-04 5.43658186e-04 2.55414133e-04 2.57925230e-04
50001815e-04 2.64633710e-04 3.71728892e-04 5.23219628e-04
08927700e-04 2.78169248e-04 3.03316115e-04 3.29103640e-04
50804495e-04 2.47784628e-04 3.43024543e-04 8.04772883e-04
68602838e-04 3.63458034e-03 3.61478735e-03 4.11608446e-04
75786430e-04 2.62410102e-04 4.64802687e-04 2.63496724e-04
33857567e-04 2.77274533e-04 5.01078667e-04 3.42090208e-04
95906573e-04 3.48992531e-04 2.74019023e-04 3.82462266e-04
49499382e-04 3.60329524e-03 4.98543051e-04 4.58754170e-04
99287784e-04 2.90066791e-04 5.11005725e-04 8.67868112e-04
03929160e-04 2.85321238e-04 1.03605731e-02 6.19344809e-04
46206722e-04 2.61928905e-03 2.89126514e-04 4.04317282e-04
14789439e-04 5.41032937e-04 4.68125716e-04 2.82483142e-04
63130780e-04 7.50205131e-04 2.66739351e-04 4.18751814e-04
22603412e-04 5.51362034e-04 7.91641879e-04 4.36816848e-04
79457580e-04 4.46525265e-04 3.12461430e-03 2.76876967e-04
58483722e-04 3.53791532e-03 2.97204974e-04 2.62999419e-04
97663321e-04 5.85841186e-03 9.09442576e-04 2.83333635e-04
26047782e-04 9.12294548e-04 3.20830967e-04 3.99222484e-04
95444999e-04 1.67166712e-03 3.23984602e-04 4.27941932e-04
11748278e-04 1.31224780e-03 1.07865097e-02 2.73367535e-04
07708117e-04 3.23280459e-04 2.69735178e-04 3.60040135e-04
85055681e-03 2.86618779e-04 2.49277549e-04 4.60375511e-04
96411485e-04 6.52882093e-04 4.90617969e-04 3.95001225e-04
21764867e-04 2.92752399e-04 4.63077330e-04 5.00309549e-04
32037995e-04 3.20749075e-04 2.54231674e-04 4.00198976e-04
54940788e-04 9.70628458e-04 2.94543782e-04 3.36519756e-04
93415300e-04 2.49376829e-04 6.95317814e-04 3.61124579e-04
69950482e-04 2.91806508e-04 2.85707017e-03 2.38613707e-03
63286803e-04 2.88384114e-04 1.85190590e-03 9.93274311e-04
94307710e-04 2.91103242e-04 3.76488541e-04 7.12226233e-04
58230854e-04 3.20857506e-03 3.27931387e-04 6.53579046e-04
80913497e-03 4.47200022e-04 4.98029532e-04 2.59734059e-04
67699033e-04 2.97946489e-04 2.51796775e-04 7.55961243e-04
70450094e-04 2.80284751e-04 7.34322115e-04 3.82171751e-04
13375643e-04 6.55482693e-04 2.78098379e-04 2.96805132e-04
86610201e-04 3.28259515e-04 4.71620905e-04 1.04241011e-03
19521175e-04 2.76714681e-04 3.06755289e-04 4.62806729e-04
69218926e-04 5.44672688e-03 3.70524290e-04 3.07698834e-03
89940345e-04 2.84599116e-04 6.65044184e-04 3.12831265e-04
56531287e-04 3.72030288e-04 3.27770091e-04 3.54258528e-04
57750234e-04 3.87354535e-04 4.36528660e-03 2.96099098e-04
38024385e-01 4.32715339e-04 2.37162395e-03 4.54278443e-04
15313448e-04 1.24983495e-03 3.78456415e-04 1.24498717e-03
83661630e-04 4.62907874e-04 4.10149561e-04 3.86241023e-04
83869049e-04 3.01381616e-03 1.83772898e-03 3.77798003e-04
53370743e-04 3.07065505e-04 2.98605919e-04 2.73322810e-04
09264921e-04 3.14676919e-04 4.63606496e-04 4.30723878e-04
55919618e-04 5.29171324e-04 4.66093611e-04 2.10386753e-03
33256599e-04 5.02692880e-04 2.87112698e-04 3.74359189e-04
12368308e-04 8.57238075e-04 2.80246652e-04 3.21797874e-04
01113123e-03 5.10846722e-03]

In diesem Modell werden 934 Datensätze mit einer Cooks Distance > 4/n als mögliche Ausreißer identifiziert. Beispielsweise könnte versucht werden, den Einfluss dieser Datensätze auf das Modell durch Datentransformation (bspw. durch Verwendung des Logarithmus) zu verringern. Allerdings hilft diese Methode der Modellperformance genau so oft, wie sie ihr schadet. Daher wird hier auf eine Datentransformation verzichtet und die Datenpunkte stattdessen aus dem Trainingsset entfernt.

#drop outliers identified by Cook's distance
train_dataset2 = train_dataset.drop(train_dataset.index[out_d])
train_dataset2.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 15405 entries, 14185 to 17879
Data columns (total 12 columns):
 #   Column                Non-Null Count  Dtype   
---  ------                --------------  -----   
 0   housing_median_age    15405 non-null  int64   
 1   total_rooms           15405 non-null  int64   
 2   total_bedrooms        15405 non-null  int64   
 3   population            15405 non-null  int64   
 4   households            15405 non-null  int64   
 5   median_income         15405 non-null  float64 
 6   median_house_value    15405 non-null  float64 
 7   ocean_proximity       15405 non-null  category
 8   people_per_household  15405 non-null  float64 
 9   rooms_per_household   15405 non-null  float64 
 10  bedrooms_per_people   15405 non-null  float64 
 11  bedrooms_per_rooms    15405 non-null  float64 
dtypes: category(1), float64(6), int64(5)
memory usage: 1.4 MB

fig = sns.scatterplot(data=train_dataset2, x='median_income', y='median_house_value')

Im Vergleich zum ursprünglichen Trainingsset, welches im Kapitel EDA ausführlich analysiert wurde, fehlen hier einige X- und Y-Ausreißer. Dadurch hat sich die waagrechte Punktelinie bei 500000$ medin_house_value verkleinert.

fig = sns.displot(data=train_dataset2, x='median_house_value')

Durch die Entfernung der Ausreißer hat sich das Histogramm von median_house_value etwas weiter normalisiert.

# Fit Model without outliers
lm = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household', data=train_dataset2).fit()
lm.summary()

OLS Regression Results
Dep. Variable:	median_house_value	R-squared:	0.741
Model:	OLS	Adj. R-squared:	0.741
Method:	Least Squares	F-statistic:	8832.
Date:	Wed, 19 Jan 2022	Prob (F-statistic):	0.00
Time:	12:37:35	Log-Likelihood:	-1.8925e+05
No. Observations:	15405	AIC:	3.785e+05
Df Residuals:	15399	BIC:	3.786e+05
Df Model:	5
Covariance Type:	nonrobust

	coef	std err	t	P>\|t\|	[0.025	0.975]
Intercept	3.596e+04	4371.548	8.226	0.000	2.74e+04	4.45e+04
ocean_proximity[T.INLAND]	-6.127e+04	1092.791	-56.071	0.000	-6.34e+04	-5.91e+04
ocean_proximity[T.NEAR COAST]	6831.6716	1113.689	6.134	0.000	4648.711	9014.633
median_income	4.676e+04	362.691	128.920	0.000	4.6e+04	4.75e+04
bedrooms_per_rooms	3.536e+05	1.2e+04	29.400	0.000	3.3e+05	3.77e+05
people_per_household	-2.532e+04	576.811	-43.895	0.000	-2.64e+04	-2.42e+04

Omnibus:	901.834	Durbin-Watson:	1.985
Prob(Omnibus):	0.000	Jarque-Bera (JB):	1114.776
Skew:	0.584	Prob(JB):	8.50e-243
Kurtosis:	3.612	Cond. No.	153.

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Durch die Entfernung der Ausreißer haben sich R-squared, adjusted R-squared und F-statistics deutlich verbessert.
AIC und BIC sind auch etwas kleiner geworden und haben sich somit verbessert.
Die p-values der Features bleiben unverändert.
Auch an den p-values für Jarque-Bera und Omnibus hat sich kaum etwas verändert.

fig = sm.graphics.plot_partregress_grid(lm)
fig.tight_layout(pad=0.2)

Die Partial Regression Plots visualisieren die Relation einer unabhängigen Variablen in Verbindung mit den anderen unabhängigen Variablen im Modell in Bezug auf die abhängige Variable.
Die Plots zeigen eine negative Abhängigkeit von T.INLAND und peope_per_household auf median_house_value. Median_income und bedrooms_per_rooms scheinen eine positive Abhängigkeit auf median_house_value zu haben. Gemäß dem Diagramm hat die Variable ocean_proximity[T.NEAR COAST] keinen Einfluss auf median_income. Dies steht im Gegensatz zum p-value 0,00 dieser Variablen, nach dessen Interpretation die Variable einen Einfluss auf die abhängige Variable hat.

fig = sm.graphics.plot_regress_exog(lm, "ocean_proximity[T.NEAR COAST]")
fig.tight_layout(pad=1.0)

Auch diese Plots lassen darauf schließen, dass die Variable _ocean_proximity[T.NEAR COAST] das Modell nicht verbessert. Sowohl Partial Regression Plot als auch CCPR Plot visualisieren eine Nicht-Linearität zwischen der Variablen und median_house_value.
Die Residuals scheinen dagegen gleichmäßig verteilt zu sein.

fig = sm.graphics.plot_regress_exog(lm, "ocean_proximity[T.INLAND]")
fig.tight_layout(pad=1.0)

Die negative Korrelation ist sowohl im Partial Regression Plot als auch im CCPR Plot zu erkennen. Die Residuals sind gleichmäßig verteilt und lassen daher auf Homoskedastizität rückschließen. Das Diagramm zu Y and Fitted vs. X ist aufgrund der vielen Datenpunkte schwer zu analysieren.

fig = sm.graphics.plot_regress_exog(lm, "people_per_household")
fig.tight_layout(pad=1.0)

Auch hier ist eine negative Abhängigkeit der abhängigen zur unabhängigen Variablen zu erkennen. Linearität scheint also vorzuliegen. Allerdings ist im Residual Plot keine konstante Varianz der Fehlerstreuung zu erkennen. Mit steigendem X-Wert verringert sich die Fehlerstreuung und liegt nur noch im positiven Bereich. Mit zunehmender people_per_household wird tendenziell ein niedrigerer median_house_value mit dem Modell hervorgesagt, als tatsächlich vorliegt. Dies spiegelt sich auch im Plot zu Y and Fitted vs. X wieder. Bei dieser Variablen könnte Heteroskedastizität vorlieren. Dies gilt es mithilfe des Breusch-Pagan Lagrange Multiplier Tests genauer zu untersuchen:

name = ['Lagrange multiplier statistic', 'p-value', 'f-value', 'f p-value']
test = sm.stats.het_breuschpagan(lm.resid, lm.model.exog)
lzip(name, test)

[('Lagrange multiplier statistic', 1043.0254994601257),
 ('p-value', 2.906700494490131e-223),
 ('f-value', 223.66770899896423),
 ('f p-value', 3.329323077460106e-231)]

Da alle Werte über 0,05 liegen, kann die Nullhypothese, Homoskedastizität, akzeptiert werden. Somit liegt keine Hetereoskedastizität (Korrelation der Fehlerterms) vor.

fig = sm.graphics.plot_regress_exog(lm, "median_income")
fig.tight_layout(pad=1.0)

Für median_income ist die Linearität im Partial Regression Plot und CCPR Plot deutlich zu erkennen. Jedoch nimmt auch hier die Streuung der Residuals im höheren X-Wertebereich ab. Heteroskedastizität wurde mithilfe des Breusch-Pagan Lagrange Multiplier Tests bereits ausgeschlossen.

Für das nächste Modell werden die Ausprägungen NEAR COST und <1H OCEAN zusammengefasst:

# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR COAST': 'NEAR COAST',
    'INLAND': 'INLAND',
    '<1H OCEAN': 'NEAR COAST' }

train_dataset.ocean_proximity = train_dataset.ocean_proximity.map(remap_cat_dict).astype('category')

Model 2 - OLS Regression¶

Nachfolgend wird ein weiteres Modell mit zusätzlicher Verwendung der Variablen housing_median_age trainiert.

# Fit Model
lm2 = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household + housing_median_age', data=train_dataset).fit()
lm2.summary()

OLS Regression Results
Dep. Variable:	median_house_value	R-squared:	0.632
Model:	OLS	Adj. R-squared:	0.631
Method:	Least Squares	F-statistic:	5598.
Date:	Wed, 19 Jan 2022	Prob (F-statistic):	0.00
Time:	12:57:18	Log-Likelihood:	-2.0552e+05
No. Observations:	16339	AIC:	4.110e+05
Df Residuals:	16333	BIC:	4.111e+05
Df Model:	5
Covariance Type:	nonrobust

	coef	std err	t	P>\|t\|	[0.025	0.975]
Intercept	-5.941e+04	4251.872	-13.972	0.000	-6.77e+04	-5.11e+04
ocean_proximity[T.NEAR COAST]	6.572e+04	1340.461	49.026	0.000	6.31e+04	6.83e+04
median_income	4.432e+04	403.597	109.805	0.000	4.35e+04	4.51e+04
bedrooms_per_rooms	2.979e+05	1.28e+04	23.242	0.000	2.73e+05	3.23e+05
people_per_household	-1.439e+04	487.311	-29.524	0.000	-1.53e+04	-1.34e+04
housing_median_age	1009.9369	45.939	21.984	0.000	919.891	1099.983

Omnibus:	5235.108	Durbin-Watson:	1.968
Prob(Omnibus):	0.000	Jarque-Bera (JB):	34340.227
Skew:	1.379	Prob(JB):	0.00
Kurtosis:	9.545	Cond. No.	766.

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

R-squared und adj. R-squared haben sich im Vergleich zum ersten Modell verbessert.
Auch F-statistics hat sich durch die Entfernung der Ausprägung von ocean_proximity und Aufnahme des Features median_housing_age verbessert im Vegleich zu Modell 1 ohne Entfernung der Cook’s Distance Outliers.
Alle übrigen Kennzahlen haben sich kaum verändert.

# obtain Cook's distance 
lm_cooksd = lm2.get_influence().cooks_distance[0]

# get length of df to obtain n
n = len(train_dataset["median_income"])
print('n = ', n)

# calculate critical d
critical_d = 4/n
print('Critical Cooks distance:', critical_d)

# identification of potential outliers with leverage
out_d = lm_cooksd > critical_d

# output potential outliers with leverage
print(train_dataset.index[out_d], "\n", 
    lm_cooksd[out_d])

n =  16339
Critical Cooks distance: 0.0002448130240528796
Int64Index([15778,  6776, 10005, 11511,  7263, 15723, 14756, 15687,  4638,
            14806,
            ...
            15809,  8866,   199, 18498,    59, 10749,  8857, 13154,  5470,
            16642],
           dtype='int64', length=855) 
 [1.33657972e-03 3.30853996e-04 2.57049009e-04 3.39660405e-04
01129532e-04 3.88602976e-04 1.90478556e-02 6.27488571e-04
74476388e-04 4.45045388e-04 8.27915922e-04 1.30029194e-03
22608518e-03 3.77876614e-04 5.13746812e-04 7.71608777e-04
53086958e-04 2.60583344e-04 3.05990146e-04 5.38968511e-04
38650064e-04 2.54349021e-04 4.32750022e-04 3.44212356e-04
61317235e-04 3.67082550e-04 2.50921062e-04 2.64010646e-04
64470925e-02 4.50645263e-04 4.32630903e-04 3.73737786e-04
94441250e-04 5.17859249e-04 1.67765900e-03 2.15661985e-03
05909753e-03 2.68151181e-04 6.29731152e-04 5.30310328e-04
08429929e-03 5.19828279e-04 9.65158028e-04 4.65438401e-03
30201917e-04 3.11323002e-04 2.63976791e-04 9.17685336e-04
92775090e-04 3.04159208e-04 4.29591546e-04 3.16585850e-04
73111270e-04 3.40214397e-04 3.69651789e-04 4.12626701e-04
19689965e-03 3.24624515e-04 6.25756905e-04 3.70903135e-04
22827592e-04 3.95367319e-04 1.05559881e-03 2.65758810e-04
47592532e-04 8.69636078e-04 2.51911954e-04 4.41768820e-04
76831782e-04 5.65039902e-04 3.33787249e-04 7.45292365e-04
02617996e-03 5.39873094e-04 4.40784778e-04 4.95448437e-04
26999608e-04 2.77340040e-04 1.84706601e-03 2.91896086e-04
69739979e-04 2.17140974e-03 3.09705912e-04 2.47578772e-04
29639650e-04 7.24422015e-04 6.79907575e-04 3.33754238e-04
15815108e-04 1.02444650e-03 3.00226007e-04 8.52987194e-04
19880125e-03 9.80927157e-04 5.65469194e-04 9.45727023e-04
21371291e-03 3.47389825e-04 6.50989809e-04 4.49479516e-04
65423758e-03 3.89116048e-04 9.94555835e-04 7.00191096e-04
18971913e-03 2.85484826e-04 2.49718013e-04 2.47484819e-04
81402090e-04 5.88963038e-04 3.45832307e-04 6.95567672e-04
35200588e-04 4.26719218e-04 2.96215870e-04 2.51307265e-04
55813701e-04 3.63345137e-04 8.04986478e-04 2.78317071e-04
96319393e-04 8.63988528e-04 1.32463248e-03 2.87769237e-04
52357194e-04 3.93632745e-04 3.92028465e-04 4.05636189e-04
50348707e-04 5.69254539e-04 4.26315065e-04 4.36699075e-02
33412362e-04 5.49398317e-04 7.26014229e-04 1.33508397e-03
01538018e-04 2.47242706e-04 5.71897176e-04 4.94021524e-04
90467300e-04 3.25255208e-04 3.04205545e-04 3.20955852e-04
00539166e-04 2.65241798e-04 3.52632314e-04 4.16500906e-04
33311859e-04 4.41249366e-04 2.86018112e-04 9.48344105e-04
47757062e-04 1.36943368e-03 5.11072467e-03 4.05563886e-04
26238220e-04 2.53555665e-04 4.54670905e-04 2.61021149e-04
67450043e-04 4.35513441e-04 3.45787798e-04 6.76666771e-04
23889590e-03 7.43737556e-04 2.79617027e-04 2.58099362e-04
38076353e-03 3.30948644e-03 3.16852303e-04 3.51331295e-04
46597653e-03 7.91184137e-04 4.22485707e-03 9.02085835e-04
95820079e-04 2.73525513e-04 3.67308349e-04 1.23677595e-03
90126981e-03 3.65851977e-04 5.99874448e-04 2.58979835e-04
28871739e-03 2.63559446e-04 4.61483732e-04 2.46606867e-02
76826150e-03 3.96818944e-03 5.81099998e-04 1.47960745e-03
60684593e-04 3.97478169e-04 2.64012779e-04 3.59647903e-04
17153655e-03 2.51225904e-04 3.34967441e-04 3.21499259e-04
71827753e-04 6.23453463e-04 1.02455371e-03 2.88236000e-04
28846792e-04 4.54854733e-04 3.02734668e-04 2.56153363e-04
02535591e-03 7.74844873e-02 2.95057072e-04 2.74284063e-04
12112604e-03 3.79537171e-04 4.34335669e-04 1.46176111e-03
81276561e-03 6.18689578e-04 3.64098379e-04 2.99725974e-04
41121691e-02 8.69079731e-04 4.63531291e-04 4.38290542e-04
95782673e-04 2.94414239e-04 1.45253916e-03 4.50655418e-04
88441108e-04 3.87161378e-04 3.66006522e-04 9.21210440e-04
62048419e-04 3.31125955e-03 4.62450237e-04 6.10613510e-04
83632763e-04 3.13751494e-04 5.73204750e-04 9.66607950e-04
66662549e-04 1.86712979e-03 2.70711139e-04 3.30637939e-04
16621222e-04 4.80582597e-03 2.74162517e-04 8.27324123e-04
37910085e-03 2.99351426e-04 3.04995191e-03 4.56892388e-04
52613038e-04 5.36569442e-04 5.52969777e-04 6.79600369e-04
71211781e-04 2.55522024e-04 4.41834960e-04 4.99009844e-04
31654890e-04 2.71205205e-03 3.09116327e-04 3.83554603e-04
33114795e-04 8.27489363e-04 3.20148895e-03 4.79803828e-03
39460248e-04 4.14850296e-04 2.98175198e-04 2.74712755e-04
68523367e-04 2.96752158e-04 3.48303317e-04 4.35518144e-04
34146132e-04 3.58313083e-04 3.62306067e-03 5.94764207e-04
82957434e-03 5.11217980e-04 3.57485411e-03 1.82299160e-03
36522236e-03 2.98418005e-04 3.56942147e-04 4.66868862e-04
08039654e-03 4.00715660e-04 3.05312182e-04 6.50790027e-04
29660141e-04 5.15219091e-04 4.36037488e-03 9.54040651e-04
30418181e-03 6.11263987e-04 1.46496999e-03 6.57213027e-03
61310109e-04 1.45812009e-03 3.38598482e-04 4.34048838e-04
71587663e-04 6.17032432e-04 6.01449533e-04 2.82616984e-04
31176835e-03 2.20882394e-03 2.90881325e-04 3.70669677e-04
27466076e-03 3.82776893e-04 3.10523565e-04 2.08058578e-03
09262135e-03 2.48706180e-04 5.54920785e-04 5.74331637e-04
73627368e-04 6.16772075e-04 4.97419818e-04 6.56238933e-04
77115720e-03 4.01074541e-04 3.94210327e-04 5.49362358e-04
58439557e-04 3.61686777e-04 2.54864542e-04 3.61748320e-04
27643619e-04 1.04546293e-03 2.83205488e-04 3.29696244e-04
67081383e-04 7.56376112e-04 2.45446931e-04 3.36356162e-04
80408253e-04 4.76466196e-04 4.74342320e-04 2.65400565e-04
40893288e-03 6.85550045e-04 4.49351575e-04 1.48871279e-03
62762831e-04 4.13256732e-04 8.10369204e-04 3.23475823e-04
19468061e-04 2.99378497e-04 7.42618687e-04 2.47048695e-04
57538872e-04 2.65767072e-04 3.98523666e-04 4.16083661e-04
19607817e-03 3.06136805e-04 7.32927370e-04 1.50379308e-03
62762594e-04 4.53491898e-04 3.78522623e-04 4.31704608e-04
75004759e-03 1.25032776e-03 3.27301325e-04 3.12266393e-04
51865473e-04 4.88383039e-04 3.13137587e-04 3.88165212e-04
38210593e-04 1.70122071e-03 2.52645185e-04 3.84139581e-04
15484246e-04 4.05445951e-03 5.90037902e-04 5.01378207e-03
75757842e-04 3.74035113e-04 1.35626379e-03 4.46171104e-03
61642764e-04 4.69177748e-04 3.63259146e-04 3.00602698e-03
98867073e-04 8.40009249e-04 5.37346944e-04 4.74872303e-03
07012536e-04 4.71755072e-03 2.79907198e-04 3.72254610e-04
46859136e-04 7.01695188e-04 3.62247955e-04 2.62301041e-04
40390396e-03 3.82269562e-04 4.59987440e-03 4.79291848e-04
23994832e-04 8.71842545e-04 7.31086891e-04 2.91781864e-04
53599743e-04 3.79510317e-04 5.19125563e-04 2.00023531e-03
88562380e-04 3.87007973e-03 6.70528630e-04 7.00383138e-04
11933556e-04 3.30141785e-04 7.90048123e-04 1.51745140e-03
73859080e-04 5.51368068e-04 3.40206222e-01 4.09083815e-03
76891564e-04 4.40940057e-04 6.65542584e-04 5.19568296e-04
39848688e-03 4.80237782e-04 2.45207300e-04 5.29121233e-03
64659410e-04 1.55945230e-03 5.47967301e-04 7.98754034e-04
95172430e-04 2.07589838e-03 3.09065176e-04 3.57159833e-04
11161736e-04 1.14456899e-03 2.32668001e-03 3.03037007e-04
16294004e-04 2.17498701e-03 3.06452930e-04 2.15205450e-02
58303543e-04 6.51516359e-04 6.42360248e-04 5.02449866e-04
80919004e-04 3.89958234e-04 1.56892155e-03 3.20696952e-04
40111749e-03 3.62088526e-04 8.85107059e-03 3.09165330e-04
37383808e-04 7.56985962e-04 5.16424609e-03 7.23450017e-03
30088735e-04 2.72886104e-04 1.36296752e-03 2.53349179e-04
25696693e-04 3.72832272e-03 4.24933142e-04 4.78915688e-04
88675846e-04 3.51261179e-04 1.90341225e-03 4.30589697e-04
55641258e-04 3.28793025e-04 7.57311788e-04 2.28638966e-03
68079850e-04 4.08505951e-04 5.71515284e-04 4.47527006e-04
77598892e-04 5.20534070e-04 2.48821737e-04 6.67810160e-04
80479901e-04 4.77746527e-04 2.58099360e-04 5.27701037e-04
01523059e-03 9.41717491e-03 1.76423553e-03 4.64299667e-04
18435904e-04 4.23968342e-04 4.18389898e-03 2.93162192e-04
30443196e-03 6.57153863e-04 4.94020048e-04 2.73543728e-04
55453033e-03 3.55310572e-04 1.66307074e-03 3.00990435e-03
59473297e-04 1.14708882e-03 3.37876852e-04 2.43977525e-03
33451538e-04 4.82836805e-04 7.05398513e-04 1.03345980e-03
28729542e-03 4.88666974e-03 8.72852756e-04 2.70601299e-04
63579349e-04 7.04879953e-04 1.70529186e-03 2.98582889e-04
80295037e-04 9.45019235e-04 2.45798045e-04 6.77037524e-04
80978771e-04 3.70850286e-04 3.85126597e-04 2.60652287e-04
03591296e-04 3.83213722e-04 2.42449874e-03 2.62531118e-04
73664514e-04 5.34530670e-04 6.26206579e-04 2.10418500e-03
13529876e-03 3.99757686e-04 4.80677677e-03 2.64539355e-04
57468005e-04 8.71438026e-04 3.05030313e-04 3.38028708e-04
24752940e-04 2.88237398e-04 6.36844746e-04 6.11920957e-04
09884554e-04 3.78816348e-03 3.18898689e-03 3.90350035e-04
21977609e-03 4.34112472e-04 7.11096691e-04 3.65261080e-04
18131957e-03 5.19152658e-04 2.25585445e-03 1.31315613e-03
28507762e-04 2.60907670e-04 2.81780237e-04 4.40587385e-04
29770417e-04 6.18218125e-04 5.46538127e-04 4.73675710e-04
33895448e-03 2.96708169e-04 5.68386272e-04 2.79242822e-04
32904360e-04 9.40540213e-04 8.73715661e-04 1.08553899e-03
69807333e-04 4.68858814e-04 2.57848006e-04 1.04619269e-03
96015914e-04 5.94117120e-04 2.48185127e-04 2.85245992e-04
28402505e-04 1.18616675e-03 9.29303308e+00 5.49886390e-04
53443430e-04 8.98135826e-04 8.49355748e-04 3.41861254e-04
02769323e-04 2.59610601e-04 6.57754109e-04 3.04390099e-04
65636458e-04 1.00797425e-03 2.89614985e-04 7.52516897e-04
52428481e-04 2.57021878e-04 3.45776714e-04 4.25440996e-04
95424860e-04 4.93146020e-04 6.31607100e-04 3.37135500e-04
85598420e-04 3.41713065e-04 3.47394867e-03 4.75862673e-04
99887261e-04 5.28215197e-04 1.94112686e+00 2.43637158e-03
97398536e-04 6.39634071e-04 2.58069091e-04 3.64716610e-04
04666472e-04 3.11696796e-04 4.26392216e-04 2.86567595e-04
63680529e-04 3.52411077e-04 3.13099501e-04 2.79049787e-04
02202980e-04 2.49901590e-04 3.82921309e-03 3.37371691e-04
91423632e-04 5.39941052e-03 4.48885326e-04 4.07443246e-04
33543785e-04 4.67380034e-04 5.37602071e-03 4.11588986e-04
19731351e-04 4.26044123e-03 5.08720585e-04 1.05689846e-03
72293013e-04 4.70739582e-04 6.20712803e-04 2.54363847e-04
82760608e-04 1.71267508e-03 6.26260569e-04 5.61326092e-04
15994858e-04 7.27738337e-04 2.63183174e-04 4.69703159e-04
54973559e-04 2.47869166e-04 2.98776045e-04 3.37824715e-04
52307238e-04 2.99002790e-04 8.14543684e-04 1.14853992e-03
06111547e-04 2.70671448e-04 3.58912851e-03 2.95878872e-03
70420200e-04 3.63572997e-04 5.58891976e-04 3.11332134e-04
47067682e-04 2.66852940e-04 2.49295725e-04 4.56368582e-04
03073713e-04 7.79036651e-04 4.45589887e-04 2.63499974e-04
71455846e-04 3.58633169e-04 6.22896251e-04 3.41261030e-03
59483680e-04 2.79288693e-04 3.28092543e-04 4.93296139e-04
18207619e-04 3.58118793e-04 9.22919457e-03 1.01698828e-03
46036129e-03 2.69571827e-04 2.69413324e-04 6.09778570e-04
16621065e-04 3.27362048e-04 7.84633461e-04 5.69964234e-04
94273342e-04 7.88776081e-04 5.73427166e-04 7.69329874e-04
63132132e-04 3.82785182e-04 4.80224041e-04 4.67819955e-03
93183311e-04 2.85653497e-04 5.26704673e-04 4.49319583e-03
60544020e-04 3.03403553e-04 3.83932886e-04 2.73482078e-04
36864469e-04 6.29663818e-03 1.93829988e-03 2.65873474e-04
29386083e-04 7.98255073e-04 3.45317099e-04 2.92108251e-04
85901435e-04 2.00085431e-03 2.63800696e-04 3.02523732e-04
95086387e-04 9.88739159e-04 9.79618143e-03 2.65917969e-04
95889252e-04 3.34933476e-04 3.33607404e-04 4.39227563e-03
37565380e-04 2.84080493e-04 2.57103220e-04 4.87948003e-04
74950736e-04 7.35561526e-04 3.28153321e-04 3.55074314e-04
34911661e-04 4.86693557e-04 3.94888958e-04 4.63609727e-04
39861959e-04 3.25750411e-04 3.37174171e-04 3.59227491e-04
73408109e-04 1.00420394e-03 3.34189302e-04 4.73270827e-04
01370224e-04 3.21689031e-04 4.10616719e-03 2.63910823e-03
82457950e-04 1.13528053e-03 1.02591202e-03 2.50920965e-04
49476617e-04 3.74311851e-04 5.73641592e-04 4.62242764e-04
82551814e-03 6.97545209e-04 4.25620530e-03 4.76380839e-04
60194360e-04 2.68900001e-04 4.65780634e-04 7.61949117e-04
74783284e-04 2.93135921e-04 3.79753843e-04 7.05919878e-04
92907105e-04 2.64060448e-04 5.82742792e-04 3.13448748e-04
92214076e-04 6.46059214e-04 4.25092730e-04 4.81368983e-04
61617768e-04 3.74740040e-04 2.69641578e-04 3.59004818e-04
92355790e-04 3.71278909e-04 6.40900628e-03 3.12586210e-04
58056383e-04 3.11797869e-03 4.47675955e-04 7.82446476e-04
59928280e-04 2.50198218e-04 5.57714432e-04 3.13549223e-04
45425403e-04 3.04159822e-04 6.08585190e-04 4.45880833e-04
86878759e-03 6.92692489e-04 9.23881597e-01 5.13395776e-04
02779635e-03 3.41076341e-04 5.53810879e-04 9.13558594e-04
36329717e-04 1.48721106e-03 1.06533957e-03 4.87951157e-04
78071534e-04 2.82509696e-04 2.84646006e-03 2.34202207e-03
97507204e-04 3.66535812e-04 2.70075220e-04 1.09193708e-03
67826847e-04 3.65926883e-04 2.71532543e-04 3.79689580e-04
95138559e-04 6.74076825e-04 3.81687946e-04 3.67945522e-04
43888506e-03 3.22455718e-04 3.78351328e-04 3.17660904e-04
11798201e-04 3.06942337e-04 9.25289483e-04 2.76054872e-04
44925121e-04 9.70486689e-04 5.76658032e-03]

Auch in Modell lm2 werden mögliche outliers und Punkte mit high leverage identifiziert und aus dem Datenset entfernt.

#drop outliers identified by Cook's distance
train_dataset2 = train_dataset.drop(train_dataset.index[out_d])
train_dataset2.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 15484 entries, 14185 to 17879
Data columns (total 12 columns):
 #   Column                Non-Null Count  Dtype   
---  ------                --------------  -----   
 0   housing_median_age    15484 non-null  int64   
 1   total_rooms           15484 non-null  int64   
 2   total_bedrooms        15484 non-null  int64   
 3   population            15484 non-null  int64   
 4   households            15484 non-null  int64   
 5   median_income         15484 non-null  float64 
 6   median_house_value    15484 non-null  float64 
 7   ocean_proximity       15484 non-null  category
 8   people_per_household  15484 non-null  float64 
 9   rooms_per_household   15484 non-null  float64 
 10  bedrooms_per_people   15484 non-null  float64 
 11  bedrooms_per_rooms    15484 non-null  float64 
dtypes: category(1), float64(6), int64(5)
memory usage: 1.4 MB

# Fit Model
lm2 = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household + housing_median_age', data=train_dataset2).fit()
lm2.summary()

OLS Regression Results
Dep. Variable:	median_house_value	R-squared:	0.746
Model:	OLS	Adj. R-squared:	0.746
Method:	Least Squares	F-statistic:	9085.
Date:	Wed, 19 Jan 2022	Prob (F-statistic):	0.00
Time:	14:04:48	Log-Likelihood:	-1.9034e+05
No. Observations:	15484	AIC:	3.807e+05
Df Residuals:	15478	BIC:	3.807e+05
Df Model:	5
Covariance Type:	nonrobust

	coef	std err	t	P>\|t\|	[0.025	0.975]
Intercept	-4.951e+04	4086.781	-12.115	0.000	-5.75e+04	-4.15e+04
ocean_proximity[T.NEAR COAST]	5.763e+04	1060.730	54.330	0.000	5.56e+04	5.97e+04
median_income	4.857e+04	368.630	131.763	0.000	4.78e+04	4.93e+04
bedrooms_per_rooms	3.637e+05	1.19e+04	30.623	0.000	3.4e+05	3.87e+05
people_per_household	-2.752e+04	578.566	-47.572	0.000	-2.87e+04	-2.64e+04
housing_median_age	923.2678	36.548	25.262	0.000	851.629	994.907

Omnibus:	1031.972	Durbin-Watson:	1.988
Prob(Omnibus):	0.000	Jarque-Bera (JB):	1351.257
Skew:	0.610	Prob(JB):	3.79e-294
Kurtosis:	3.778	Cond. No.	912.

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Die Modellperformance liegt nun bei 74,6% und ist somit bisher der beste erzielte Wert der Kennzahl R-squared.
F-statistics hat sich ebenfalls verbessert.
Alle anderen Werte liegen in einem ähnlichen Wertebereich.

fig = sm.graphics.plot_partregress_grid(lm2)
fig.tight_layout(pad=0.2)

Alle Features zeigen einen Einfluss auf median_house_value.

fig = sm.graphics.plot_regress_exog(lm2, "housing_median_age")
fig.tight_layout(pad=1.0)

Das Partial Regression Plot und CCPR Plot zeigen eine geringe lineare Korrelation von housing_median_age in Verbindung mit den anderen unabhängigen Variablen auf median_house_value.
Das Residual Plot zeigt eine deutliche Homoskedastizität.

fig = sm.graphics.plot_regress_exog(lm2, "ocean_proximity[T.NEAR COAST]")
fig.tight_layout(pad=1.0)

Durch die Zusammenlegung der Ausprägungen <1H OCEAN und NEAR COST hat ds Feature nun eine positive lineare Abhängigkeit zu median_house_value. Die Fehlerstreuung ist bei Distrikten in Meerehöhe etwas größer. Auch hier wird vorsichtshalber der Breusch-Pagan Lagrange Multiplier Test zum Ausschluss von Heteroskedastizität durchgeführt.

name = ['Lagrange multiplier statistic', 'p-value', 'f-value', 'f p-value']
test = sm.stats.het_breuschpagan(lm2.resid, lm2.model.exog)
lzip(name, test)

[('Lagrange multiplier statistic', 1095.2984925363032),
 ('p-value', 1.3939938205969558e-234),
 ('f-value', 235.64364107050312),
 ('f p-value', 2.5613416182330695e-243)]

Da alle p-values über 0,05 liegen, kann Heteroskedastizität ausgeschlossen werden.

In der linearen Regression muss eine Colinearität der unabhängigen Variablen unbedingt vermieden werden. Mithilfe des Variance Inflation Factors (VIF) kann der Grad der Colinearität der im Modell verwendeten Variablen berechnet werden. Der kleinstmögliche Wert ist 1. Problematisch wird es ab Werten > 5.

# choose features and add constant
features = add_constant(train_dataset2[['median_income', 'bedrooms_per_rooms', 'people_per_household', 'housing_median_age']])
# create empty DataFrame
vif = pd.DataFrame()
# calculate vif
vif["VIF Factor"] = [variance_inflation_factor(features.values, i) for i in range(features.shape[1])]
# add feature names
vif["Feature"] = features.columns

vif.round(2)

C:\ProgramData\Anaconda3\lib\site-packages\statsmodels\tsa\tsatools.py:142: FutureWarning:

In a future version of pandas all arguments of concat except for the argument 'objs' will be keyword-only

	VIF Factor	Feature
0	83.61	const
1	1.80	median_income
2	1.79	bedrooms_per_rooms
3	1.01	people_per_household
4	1.03	housing_median_age

Der Wert der Konstanten kann ignoriert werden. Ansonsten ist die Colinearität der Variablen eher gering.

Fazit zu Modellen 1 und 2:
Die Modelle lm und lm2 sind sehr ähnlich. R-squared ist bei lm2 ca. 1% besser, auch der Wert von F-statistics hat sich von lm zu lm2 verbessert. AIC und BIC haben sich dagegen minimal verschlechtert Das Partial Regression Plot der Variablen median_housing_age zeigt zwar nur eine geringe Linearität, dennoch scheint die Aufnahme der Variablen das Modell etwas zu verbessern. Daher verwenden die nachfolgenden Modelle die in lm2 verwendeten Features.

Model 2 - Lasso Regression¶

Nachfolgend wird eine Lasso Regression mit statsmodels fit_regularized() Methode trainiert. Der Wert von alpha wird auf 5 gesetzt.

lasso1 = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household + housing_median_age', data=train_dataset2).fit_regularized(method='sqrt_lasso', alpha=5.0, L1_wt=1.0, max_iter=1000)

print("Parameters: ", lasso1.params)

Parameters:  Intercept                        -44272.058534
ocean_proximity[T.NEAR COAST]     58127.019057
median_income                     48170.671059
bedrooms_per_rooms               346772.197084
people_per_household             -27621.154742
housing_median_age                  917.516632
dtype: float64

Da die Summary() Methode für statsmodels fit_regularized() Medelle nicht verfügbar ist, müssen wichtige Kennzahlen selbst berechnet werden.

# calculate the mean 
train_dataset2['average'] = train_dataset2["median_house_value"].mean()

# calculate error (observation - average) and assign it to dataframe
train_dataset2 = train_dataset2.assign(error = (train_dataset2['median_house_value'] - train_dataset2['average']))

# calculate squared error and assign it to dataframe
train_dataset2 = train_dataset2.assign(error_sq = (train_dataset2['median_house_value'] - train_dataset2['average'])**2)


# obtain predicted value from model lasso
train_dataset2['pred'] = lasso1.predict()

# obtain the residuals from statsmodel (resid)
train_dataset2['error_2'] = train_dataset2['median_house_value'] - train_dataset2['pred']

# square the residuals 
train_dataset2['error_sq_2'] = train_dataset2['error_2']**2

# show df
train_dataset2.head(5)

	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	people_per_household	rooms_per_household	bedrooms_per_people	bedrooms_per_rooms	average	error	error_sq	pred	error_2	error_sq_2
14185	37	2176	418	1301	375	2.8750	98900.0	NEAR COAST	3.469333	5.802667	0.321291	0.192096	196573.25704	-97673.25704	9.540065e+09	157080.171533	-58180.171533	3.384932e+09
6125	20	3158	684	2396	713	3.5250	153000.0	NEAR COAST	3.360449	4.429173	0.285476	0.216593	196573.25704	-43573.25704	1.898629e+09	184295.786398	-31295.786398	9.794262e+08
14095	11	2393	726	1905	711	1.3448	91300.0	NEAR COAST	2.679325	3.365682	0.381102	0.303385	196573.25704	-105273.25704	1.108246e+10	119926.953213	-28626.953213	8.195025e+08
14359	52	1260	202	555	209	7.2758	345200.0	NEAR COAST	2.655502	6.028708	0.363964	0.160317	196573.25704	148626.74296	2.208991e+10	394291.589328	-49091.589328	2.409984e+09
18004	32	2930	481	1336	481	6.4631	344100.0	NEAR COAST	2.777547	6.091476	0.360030	0.164164	196573.25704	147526.74296	2.176414e+10	334755.756942	9344.243058	8.731488e+07

# calculate sum of squared error (which is in case of the mean the total error)
TSS = train_dataset2.error_sq.sum()
# print output
print('Sum of squared error (TSS) of model 1:', TSS)

# Sum of squared residuals (SS_R)
SSR = train_dataset2['error_sq_2'].sum()
print('Sum of squared residuals (SSR) of model 1:',SSR)

Sum of squared error (TSS) of model 1: 169695617828075.0
Sum of squared residuals (SSR) of model 1: 43134248561958.78

# Explained sum of squares  (SS_M = TSS - SS_R)
SSM = TSS - SSR
print('SSM of Lasso Model: ', SSM)

# R_Squared: explained sum of squared residuals
R_squared = SSM / TSS
print('R squared Lasso Model:', R_squared)

SSM of Lasso Model:  126561369266116.22
R squared Lasso Model: 0.7458140103201739

R-squared ist durch die Regularisierung nochmal etwas besser geworden.
Für einen optimalen Wert von alpha müsste eigentlich eine Cross-Validation durchgeführt werden. Darauf wird an dieser Stelle jedoch verzichtet.

Model 3 - Cubic Spline¶

Nun wird eine Cubic Spline unter Verwendung von Patsy und Statsmodels trainiert. Eine Verwendung von mehreren unabhängigen Variablen ist dabei nicht ohne Weiteres möglich. Daher wird nachfolgend lediglich die Variable mit der höchsten Korrelation zu median_house_value verwendet. Dabei handelt es sich um median_income.

X_train = train_dataset2['median_income']
y_train = train_dataset2['median_house_value']
X_test = test_dataset['median_income']
y_test = test_dataset['median_house_value']

# Generating cubic spline with 3 knots at 25, 40 and 60
transformed_x = dmatrix(
            "bs(train, knots=(3,6,9), degree=3, include_intercept=False)", 
                {"train": X_train},return_type='dataframe')

# Fitting generalised linear model on transformed dataset
spline1 = sm.GLM(y_train, transformed_x).fit()

spline1.predict()
spline1.summary()

Generalized Linear Model Regression Results
Dep. Variable:	median_house_value	No. Observations:	15484
Model:	GLM	Df Residuals:	15477
Model Family:	Gaussian	Df Model:	6
Link Function:	identity	Scale:	4.5931e+09
Method:	IRLS	Log-Likelihood:	-1.9421e+05
Date:	Wed, 19 Jan 2022	Deviance:	7.1087e+13
Time:	14:14:30	Pearson chi2:	7.11e+13
No. Iterations:	3
Covariance Type:	nonrobust

	coef	std err	z	P>\|z\|	[0.025	0.975]
Intercept	1.08e+05	8224.837	13.135	0.000	9.19e+04	1.24e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[0]	-3.733e+04	1.09e+04	-3.439	0.001	-5.86e+04	-1.61e+04
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[1]	7.416e+04	7387.742	10.038	0.000	5.97e+04	8.86e+04
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[2]	1.573e+05	9964.123	15.782	0.000	1.38e+05	1.77e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[3]	4.465e+05	1.06e+04	42.023	0.000	4.26e+05	4.67e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[4]	3.512e+05	1.97e+04	17.846	0.000	3.13e+05	3.9e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[5]	4.248e+05	4.12e+04	10.307	0.000	3.44e+05	5.06e+05

# obtain predicted value from model spline
train_dataset2['pred_spline'] = spline1.predict()

# obtain the residuals from statsmodel (resid)
train_dataset2['error_2_spline'] = train_dataset2['median_house_value'] - train_dataset2['pred_spline']

# square the residuals 
train_dataset2['error_sq_2_spline'] = train_dataset2['error_2_spline']**2

	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	people_per_household	rooms_per_household	...	bedrooms_per_rooms	average	error	error_sq	pred	error_2	error_sq_2	pred_spline	error_2_spline	error_sq_2_spline
14185	37	2176	418	1301	375	2.8750	98900.0	NEAR COAST	3.469333	5.802667	...	0.192096	196573.25704	-97673.25704	9.540065e+09	157080.171533	-58180.171533	3.384932e+09	154242.775962	-55342.775962	3.062823e+09
6125	20	3158	684	2396	713	3.5250	153000.0	NEAR COAST	3.360449	4.429173	...	0.216593	196573.25704	-43573.25704	1.898629e+09	184295.786398	-31295.786398	9.794262e+08	183202.797434	-30202.797434	9.122090e+08
14095	11	2393	726	1905	711	1.3448	91300.0	NEAR COAST	2.679325	3.365682	...	0.303385	196573.25704	-105273.25704	1.108246e+10	119926.953213	-28626.953213	8.195025e+08	96483.272131	-5183.272131	2.686631e+07
14359	52	1260	202	555	209	7.2758	345200.0	NEAR COAST	2.655502	6.028708	...	0.160317	196573.25704	148626.74296	2.208991e+10	394291.589328	-49091.589328	2.409984e+09	372203.127060	-27003.127060	7.291689e+08
18004	32	2930	481	1336	481	6.4631	344100.0	NEAR COAST	2.777547	6.091476	...	0.164164	196573.25704	147526.74296	2.176414e+10	334755.756942	9344.243058	8.731488e+07	319584.526302	24515.473698	6.010085e+08

5 rows × 21 columns

# Sum of squared residuals (SSR)
SSR_spline = train_dataset2['error_sq_2_spline'].sum()
print('Sum of squared residuals (SSR) of Cubic Spline:',SSR_spline)

# Explained sum of squares  (SS_M = TSS - SS_R)
SSM_spline = TSS - SSR_spline
print('SSM of Cubic Spline: ', SSM_spline)

# R_Squared: explained sum of squared residuals
R_squared_spline = SSM_spline / TSS
print('R squared of Cubic Spline:', R_squared_spline)

Sum of squared residuals (SSR) of Cubic Spline: 71087433258433.94
SSM of Cubic Spline:  98608184569641.06
R squared of Cubic Spline: 0.581088573952126

R squared ist schlechter geworden als in den vorherigen Modellen. Allerdings liegt dies darin begründet, dass nur ein Feature verwendet wurde. Ein direkter Vergleich ist also nicht möglich.

# Create observations
xp = np.linspace(X_test.min(),X_test.max(), 100)
# Make some predictions
pred = spline1.predict(dmatrix("bs(xp, knots=(3, 6, 9), include_intercept=False)", {"xp": xp}, return_type='dataframe'))

# plot
sns.scatterplot(x=X_train, y=y_train)

plt.plot(xp, pred, label='Cubic spline with degree=3 (3 knots)', color='orange')
plt.legend();

Model 4 - Natural Spline¶

transformed_x3 = dmatrix("cr(train,df = 3)", {"train": X_train}, return_type='dataframe')

spline2 = sm.GLM(y_train, transformed_x3).fit()

# obtain predicted value from model spline
train_dataset2['pred_spline'] = spline2.predict()

# obtain the residuals from statsmodel (resid)
train_dataset2['error_2_spline'] = train_dataset2['median_house_value'] - train_dataset2['pred_spline']

# square the residuals 
train_dataset2['error_sq_2_spline'] = train_dataset2['error_2_spline']**2

# Sum of squared residuals (SSR)
SSR_spline = train_dataset2['error_sq_2_spline'].sum()
print('Sum of squared residuals (SSR) of Natural Spline:',SSR_spline)

# Explained sum of squares  (SS_M = TSS - SS_R)
SSM_spline = TSS - SSR_spline
print('SSM of Natural Spline: ', SSM_spline)

# R_Squared: explained sum of squared residuals
R_squared_spline = SSM_spline / TSS
print('R squared of Natural Spline:', R_squared_spline)

Sum of squared residuals (SSR) of Natural Spline: 72042025549993.14
SSM of Natural Spline:  97653592278081.86
R squared of Natural Spline: 0.5754632531349064

R-squared ist im Vergleich zur Cubic Spline nochmal gesunken. Auch hier wird im Modell nur eine unabhängige Variable verwendet.

# Make predictions
pred = spline2.predict(dmatrix("cr(xp, df=3)", {"xp": xp}, return_type='dataframe'))

# plot
sns.scatterplot(x=X_train, y=y_train)
plt.plot(xp, pred, color='orange', label='Natural spline with df=3')
plt.legend()

<matplotlib.legend.Legend at 0x222f091b2b0>

Scikit Learn¶

Data Preprocessing Pipeline¶

Initialisierung der Preprocessing Pipeline. Mithilfe einer Preprocessing Pipelines müssen die Daten nicht manuell transformiert werden, da diese Aufgaben von der Pipeline übernommen werden.
Für numerische Variablen ist dies beispielsweise die Normalisierung der Variablen, sowie die Belegung nicht befüllter Datensätze mit dem Durchschnittswert des Features.
Bei kategorialen Variablen werden fehlnde Werte durch eine Konstante belegt und automatisiert Dummy Variablen für die verschiedenen Ausprägungen der Variablen erzeugt.

# for numeric features
numeric_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='median')),
    ('scaler', StandardScaler()) #Normalisierung (X=(X-Mittelwert) / Standardabweichung)
    ])

# for categorical features  
categorical_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='constant', fill_value='missing')), #Konstante bei fehlenden Werten reinschreiben
    ('onehot', OneHotEncoder(handle_unknown='ignore')) #pro Featureausprägung eine Spalte, zutreffendes Feature hat dann Wert 1 und die anderen Wert 0
    ])

# Pipeline
preprocessor = ColumnTransformer(transformers=[
    ('num', numeric_transformer, selector(dtype_exclude="category")),
    ('cat', categorical_transformer, selector(dtype_include="category"))
        ])

Datenvorbereitung¶

Für Scikit Learn muss die abhängige Variable von den unabhängigen getrennt werden. Als Datenbasis wird hier auf train_dataset2 zugegriffen. In diesem Datenset wurden die Outlier nach Cook’s Distance bereits entfernt. Als Features sollen die durch die Filter Methode identifizierten und bereits in statsmodels verwendeten Features genutzt werden.

# create label y
y = train_dataset2['median_house_value']

# create features
X= train_dataset2[['median_income', 'ocean_proximity', 'bedrooms_per_rooms', 'people_per_household', 'housing_median_age']]

Scikit Learn bietet eine eigene Methode zum Datensplitting. Das Datenset train_dataset2 aus statsmodels wird nun nochmal in Test- und Trainingsdaten aufgeteilt. Dadurch ist eine Validierung der Modelle bereits mit diesem Testset möglich, bevor am Ende dieses Notebooks mit den eigentlichen Testdaten das beste Modell evaluiert wird.

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=10)

print("Training size:", len(X_train))
print("Test size:", len(X_test))

Training size: 12387
Test size: 3097

Model 5 - Simple Regression¶

Nun wird ein Pipeline Object erzeugt, welches zunächst die Data Preprocessing Pipeline ausführt und anschließend das Regressionsmodell.

# Create pipeline with model
lm_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lm', LinearRegression())
                        ])

# Fit model
lm_pipe.fit(X_train, y_train)

Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F09495B0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F0949940>)])),
                ('lm', LinearRegression())])

In der Variable features_names werden die Namen der Features gespeichert.

feature_names = np.concatenate((X_train.drop(['ocean_proximity'], axis=1).columns.to_numpy(), lm_pipe.named_steps['preprocessor'].transformers_[1][1]['onehot'].get_feature_names_out()))
feature_names

array(['median_income', 'bedrooms_per_rooms', 'people_per_household',
       'housing_median_age', 'x0_INLAND', 'x0_NEAR COAST'], dtype=object)

Der eben erzeugte feature_names Array kann nun mit den Koeffizienten des Modells ergänzt werden, um die Koeffizienten je Feature auszugeben.

# Obtain model coefficients and names
print('Koeffizienten und Namen: \n', list(zip(lm_pipe.named_steps['lm'].coef_, feature_names)))

Koeffizienten und Namen: 
 [(84852.55471482786, 'median_income'), (18735.049763174888, 'bedrooms_per_rooms'), (-19963.944594294095, 'people_per_household'), (11488.308364831402, 'housing_median_age'), (-28587.94998804061, 'x0_INLAND'), (28587.949988040604, 'x0_NEAR COAST')]

Das trainierte Modell wird nun für die Berechnung der median_house_values der Trainingsdaten verwendet.

y_pred = lm_pipe.predict(X_train)
y_pred

array([113004.32375212, 163416.09509259, 104877.65553309, ...,
        45560.80974401, 306094.8367317 , 117268.25027873])

# Metrics trainings set
print('R squared training set:', round(skm.r2_score(y_train, y_pred)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, y_pred), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, y_pred, squared=False), 2))

R squared training set: 74.47
MSE training set: 2777618420.95
RMSE training set: 52703.12

Das Modell hat einen ähnlichen R-squared wie die OLS Regression in statsmodels. Der Root Mean Squared Error des Modells beträgt 52703,12$. Das bedeutet, dass der eigentliche Wert von median_house_value im Schnitt um diesen Betrag von dem durch das Modell berechneten Ergebnis abweicht.

Im Folgenden wird das Modell mit den im Scikit Learn Datensplit erzeugten Validierungsdaten getestet.

# Validation data
pred = lm_pipe.predict(X_test)
r2_test = round(skm.r2_score(y_test, pred)*100, 2)
print('R squared validation set:', r2_test)
mse_test =skm.mean_squared_error(y_test, pred)
print('MSE validation set', round(mse_test, 2))
print('RMSE training set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))

R squared validation set: 75.04
MSE validation set 2816769494.14
RMSE training set: 53073.25

Der R-squared Wert hat sich mit den Validierungsdaten verbessert. Dies ist sehr ungewöhnliche, da sich unter Verwendung unbekannter Variablen die Kennzahlen normalerweise etwas verschlechtern. Allerdings wurde mit diesen Validierungsdaten bereits EDA, Feature Engineering und auch Cook’s Distance Outlier Entfernung durchgeführt. Die Daten sind somit nicht wirklich unbekannt.
Der durchschnittliche Fehler MSE und auch die durchschnittliche Abweichung von berechneten zu tatsächlichen Modellen sind etwas schlechter geworden. Allerdings nur geringfügig.

Model 6 - Lasso¶

vorgegebenes Alpha¶

In diesem Unterkapitel wird eine Lasso Regression für das Modell durchgeführt. Der Wert für alpha wird initial mit 1 belegt.

# Create pipeline with model
lasso_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lasso', Lasso(alpha=1, max_iter=10000))
                        ])

lasso_pipe.fit(X_train, y_train)

Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F09495B0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F0949940>)])),
                ('lasso', Lasso(alpha=1, max_iter=10000))])

# Metrics trainings set
pred_train = lasso_pipe.predict(X_train)
print('R squared training set:', round(skm.r2_score(y_train, pred_train)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, pred_train), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, pred_train, squared=False), 2))

R squared training set: 74.47
MSE training set: 2777618427.33
RMSE training set: 52703.12

R-squared der Lasso Regression hat sich im Vergleich zur Simple Regression im Modell zuvor nicht verändert. Auch MSE und RMSE sind gleich geblieben. Dies ist auf den initial gewählten, sehr kleinen Wert für alpha von 1 zurückzuführen. Ein geringer Wert in Alpha verändert die Höhe der Koeffizienten im Modell kaum. Daher bleiben auch die Berechnungen des Modells unverändert und somit die Metriken.

# Metrics validation set
pred = lasso_pipe.predict(X_test)
print('R squared validation set:', round(skm.r2_score(y_test, pred)*100, 2))
print('MSE validation set:', round(skm.mean_squared_error(y_test, pred), 2))
print('RMSE validation set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))

R squared validation set: 74.67
MSE validation set: 2877369775.47
RMSE validation set: 53641.12

Auch die Metriken des Validation Set haben sich nicht verändert im Vergleich zur Simple Regression.

# Obtain model coefficients and names
print('Koeffizienten und Namen: \n', list(zip(lasso_pipe.named_steps['lasso'].coef_, feature_names)))

Koeffizienten und Namen: 
 [(84852.14248638194, 'median_income'), (18734.010558726382, 'bedrooms_per_rooms'), (-19962.92053743495, 'people_per_household'), (11487.952227452783, 'housing_median_age'), (-57171.981131734494, 'x0_INLAND'), (3.3287087880073534e-11, 'x0_NEAR COAST')]

# get absolute values of coefficients
importance = np.abs(lasso_pipe.named_steps['lasso'].coef_)

sns.barplot(x=importance, 
            y=feature_names);

Das Diagramm visualiert die Feature Importance. Je länger der Balken, desto eine höhere Gewichtung bekommt das Feature im Modell. Die Feature Importance ergibt sich auch den Koeffizienten der Features. median_income ist das wichtigste Feature im Modell, wie an dem großen Balken sowie dem Koeffizienten des Features zu erkennen.
Die Ausprägung NEAR COAST der Variablen ocean_proximity wurde von der Lasso Regression aus der Regression fast entfernt. Dieses Verhalten ist bei Lasso gewünscht. Durch die L1-Regularisierung mit alpha führt der Algorithmus automatisch eine Feature Selection durch. Dabei werden relevante Features durch höhere Koeffizienten stärker gewichtet und unwichtige Features durch kleinere Koeffizienten schwächer gewichtet bzw. sogar komplett aus dem Modell entfernt, indem der Koeffizient auf den Wert 0 reduziert wird. Durch diese Methode soll ein overfitting des Modells durch zu viele Variablen verhindert werden.

k-fold cross validation zur Ermittlung von alpha¶

from sklearn.linear_model import LassoCV

# Lasso with 5 fold cross-validation
model_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lassoCV', LassoCV(cv=5, random_state=0, max_iter=10000))
                        ])
# Fit model
model_pipe.fit(X_train, y_train)

Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F09495B0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F0949940>)])),
                ('lassoCV', LassoCV(cv=5, max_iter=10000, random_state=0))])

#Show best value of penalization chosen by cross validation:
alpha = model_pipe.named_steps['lassoCV'].alpha_
print('Best alpha chosen by cross validation:', alpha)

Best alpha chosen by cross validation: 79.0619710978722

Auf Basis von k-fold cross validation wurde der beste Wert für Alpha ermittelt. Dieser liegt bei 79,06. Bei k-fold cross validation wird das Trainingsdatenset nochmal in k kleinere Datensets unterteilt. Eines dieser kleineren Datensets wird zur Seite gelegt und das Modell mithilfe der anderen Datensätzen und einem Wert für Alpha trainiert. Anschließend wird mit dem zur Seite gelegten kleinen Datenset eine Validierung des Modells mit den unbekannten Daten durchgeführt.
Dieser Vorgang wird für jedes der k - Datensets für verschiedene Werte von Alpha wiederholt. Das Alpha mit dem besten Ergebnis bei der Validierung wird hier ausgegeben.

Mit dem zuvor ermittelten, besten Wert von Alpha wird nun noch ein finales Lasso Modell auf Basis aller Trainingsdaten trainiert.

# Set best alpha
lasso_best_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lasso_best', Lasso(alpha=alpha))
                        ])
lasso_best_pipe.fit(X_train, y_train)
lasso_best = lasso_best_pipe.named_steps['lasso_best']

print('Model coefficients and names: \n', list(zip(lasso_best.coef_, feature_names)))

Model coefficients and names: 
 [(84682.01192725677, 'median_income'), (18560.71403670962, 'bedrooms_per_rooms'), (-19890.4973463913, 'people_per_household'), (11429.85201175577, 'housing_median_age'), (-56979.210567451206, 'x0_INLAND'), (5.031471360334191e-11, 'x0_NEAR COAST')]

Durch den größeren Wert von Alpha wurden die Koeffizienten der einzelnen Features etwas verkleinert. median_income hatte zuvor bei alpha=1 beispielsweise einen Wert von 84852,14. Nun liegt der Wert bei 84682,01.

# get absolute values of coefficients
importance = np.abs(lasso_best.coef_)

sns.barplot(x=importance, 
            y=feature_names);

An der Reihenfolge der Feature Importance hat sich durch die Veränderung von Alpha wenig verändert.

# Metrics trainings set
pred_train = lasso_best_pipe.predict(X_train)
print('R squared training set:', round(skm.r2_score(y_train, pred_train)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, pred_train), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, pred_train, squared=False), 2))

R squared training set: 74.47
MSE training set: 2777671055.36
RMSE training set: 52703.62

# Metrics tests set
pred = lasso_best_pipe.predict(X_test)
print('R squared validation set:', round(skm.r2_score(y_test, pred)*100, 2))
print('MSE validation set:', round(skm.mean_squared_error(y_test, pred), 2))
print('RMSE validation set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))

R squared validation set: 75.03
MSE validation set: 2817208350.85
RMSE validation set: 53077.38

Da sich die Koeffizienten durch den neuen Wert von Alpha kaum verändert haben, sind auch die Kennzahlen zur Performance Messung des Modells unverändert.

Model 7 - Natural Splines¶

In diesem Modell wird eine scikit learn spline in Kombination mit einer Data Preprocessing Pipeline trainiert.

# use a spline wit 4 knots and 3 degrees with a ridge regressions
spline = make_pipeline(SplineTransformer(n_knots=4, degree=3), 
                       Ridge(alpha=1))

# Integrate Preprocessor Pipeline
spline_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('spline', spline)
                        ])
                     
spline_pipe.fit(X_train, y_train)

Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F09495B0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F0949940>)])),
                ('spline',
                 Pipeline(steps=[('splinetransformer',
                                  SplineTransformer(n_knots=4)),
                                 ('ridge', Ridge(alpha=1))]))])

# Metrics trainings set
pred_train = spline_pipe.predict(X_train)
print('R squared training set:', round(skm.r2_score(y_train, pred_train)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, pred_train), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, pred_train, squared=False), 2))

R squared training set: 75.36
MSE training set: 2680356009.84
RMSE training set: 51772.15

# Metrics tests set
pred = spline_pipe.predict(X_test)
print('R squared validation set:', round(skm.r2_score(y_test, pred)*100, 2))
print('MSE validation set:', round(skm.mean_squared_error(y_test, pred), 2))
print('RMSE validation set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))

R squared validation set: 76.19
MSE validation set: 2687000792.6
RMSE validation set: 51836.29

Durch die Verwendung der Spline hat sich das Modell in Bezug auf R-sqared und auch RMSE nochmal verbessert. Das Modell deckt nun 75,36% der Varianz (Testdaten) bzw. 76,19% der Validationsdaten ab. Die durchschnittliche Abweichung der tatsächlichen Werten zu den berechneten Werten beträgt ca. 51800$.

# Create observations
xp = pd.DataFrame(data=np.linspace(X_test.min(),X_test.max(), 100), columns = X_train.drop(['ocean_proximity'], axis=1).columns)
xp = xp.assign(ocean_proximity=lambda xp: 'INLAND')


# Make some predictions
pred_new = spline_pipe.predict(xp)

# plot
sns.scatterplot(x=X_train['median_income'], y=y_train)

plt.plot(xp['median_income'], pred_new, label='Cubic spline with degree=3', color='orange')
plt.legend();

C:\Users\Admin\AppData\Local\Temp/ipykernel_5788/3121680389.py:2: FutureWarning:

Dropping of nuisance columns in DataFrame reductions (with 'numeric_only=None') is deprecated; in a future version this will raise TypeError.  Select only valid columns before calling the reduction.

Evaluierung der besten 2 Modelle mit Test Datenset¶

In diesem Kapitel wird eine Evaluierung der zwei besten Modellen mithilfe der zu Beginn beiseite gelegten Testdaten durchgeführt.
Zunächst müssen die Testdaten transformiert werden, sodass sie in gleicher Form wie die Trainingsdaten vorliegen und von den Modellen verwendet werden können.

# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR OCEAN': 'NEAR COAST',
    'NEAR BAY': 'NEAR COAST',
    'INLAND': 'INLAND',
    '<1H OCEAN': 'NEAR COAST',
    'ISLAND' : 'NEAR COAST'}

test_dataset.ocean_proximity = test_dataset.ocean_proximity.map(remap_cat_dict).astype('category')

y_test = test_dataset['median_house_value']

X_test = test_dataset[['median_income', 'ocean_proximity', 'bedrooms_per_rooms', 'people_per_household', 'housing_median_age']]

Evaluation Statsmodels OLS lm2¶

# Metrics tests set
X_test = X_test.assign(pred_test1=lambda X_test: lm2.predict(X_test))
r2_test1 = round(skm.r2_score(y_test, X_test.pred_test1)*100, 2)
print('R squared Statsmodels OLS with Test Set:', r2_test1)
print('MSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test, X_test.pred_test1), 2))
print('RMSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test, X_test.pred_test1, squared=False), 2))

R squared Statsmodels OLS with Test Set: -10.35
MSE Statsmodels OLS with Test Set: 14459492747.15
RMSE Statsmodels OLS with Test Set: 120247.63

R-squared ist in Model lm2 negativ. Somit generalisiert das Modell schlecht auf unbekannte Testdaten und stellt sogar eine Verschlechterung im Vergleich zur Verwendung des Durchschnitts als Modell dar.
Grund dafür könnte ein Overfitting während dem Modellierungsprozess sein.
Root Mean Squared Error und somit eine durchschnittliche Abweichung von tatsächlichem und berechnetem median_house_value in Höhe von 120247,63,95$ ist ebenso kein zufriedenstellendes Ergebnis.

Die berechneten Werte sollen in einem Scatterplot den tatsächlichen Werten gegenübergestellt werden. Aufgrund der Datenmenge sollen nur einige Datensätze visualisiert werden. Mithilfe eines Data Splittings werden zufällig einige Datenpunkte ausgewählt.

X_notplot, X_plot, y_notplot, y_plot = train_test_split(X_test, y_test, test_size=0.1, random_state=10)

# plot
sns.scatterplot(x=X_plot['median_income'], y=y_plot, label='actual')

plt.plot(X_plot['median_income'], X_plot['pred_test1'], label='predicted', color='orange', marker='.', linestyle='')
plt.legend();

Im Plot ist zu erkennen, dass bei bei einem median_income von einen Datenpunkt gibt, dessen median_house_value negativ geschätzt wurde. Dieser Wert ist unplausibel, da ein Haus keinen negativen Wert hat.
Den Koeffizienten des Modells lm2 ist zu entnehmen, dass nur das Feature people_per_household eine negative Abhängigkeit auf median_house_value hat. Die Ursache für einen negativen Wert wird also durch einen hohen Wert im Feature people_per_household verursacht.

lm2.params

Intercept                        -49511.050850
ocean_proximity[T.NEAR COAST]     57629.618091
median_income                     48571.655721
bedrooms_per_rooms               363707.323412
people_per_household             -27523.527384
housing_median_age                  923.267802
dtype: float64

test_dataset.loc[test_dataset['people_per_household'] > 100]

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	people_per_household	rooms_per_household	bedrooms_per_people	bedrooms_per_rooms
13034	-121.15	38.69	52	240	44	6675	29	6.1359	225000.0	INLAND	230.172414	8.275862	0.006592	0.183333

Tatsächlich existiert ein Datensatz mit über 100 people_per_household. Dieser ist vermutlich für den negativen R-squared Wert verantwortlich. Daher wird dieser testweise entfernt und R-squared erneut berechnet.

X_test2 = X_test.drop(13034)
y_test2 = y_test.drop(13034)

# Metrics tests set
X_test2 = X_test2.assign(pred_test1=lambda X_test2: lm2.predict(X_test2))
r2_test1 = round(skm.r2_score(y_test2, X_test2.pred_test1)*100, 2)
print('R squared Statsmodels OLS with Test Set:', r2_test1)
print('MSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test2, X_test2.pred_test1), 2))
print('RMSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test2, X_test2.pred_test1, squared=False), 2))

R squared Statsmodels OLS with Test Set: 61.36
MSE Statsmodels OLS with Test Set: 5064559536.22
RMSE Statsmodels OLS with Test Set: 71165.72

Wie vermutet, lag der Grund für den negativen Wert in diesem Ausreiser. Jedoch darf auf den Testdaten gemäß der Literatur keine EDA oder Feature Engineering betrieben werden.
Das Modell lm2 ist daher kein gutes Modell, welches ausreichend unbekannte Daten generalisiert.

Evaluation Scikit Learn Spline¶

#Calculate R2
X_test = X_test.assign(pred_test2=lambda X_test: spline_pipe.predict(X_test))
r2_test = round(skm.r2_score(y_test, X_test.pred_test2)*100, 2)
print('R squared test set:', r2_test)
mse_test =skm.mean_squared_error(y_test, X_test.pred_test2)
print('MSE test set', round(mse_test, 2))
print('RMSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test, X_test.pred_test2, squared=False), 2))

R squared test set: 65.54
MSE test set 4515256968.94
RMSE Statsmodels OLS with Test Set: 67195.66

Dieses Modell berechnet aus dem kompletten Test Datenset einen R-squared Wert von 65,54%. Damit wird das zu Beginn gesetzte Ziel, einen R-squared von 70% zu erreichen unter Verwendung der Testdaten leider verfehlt. Unter Verwendung der Trainingsdaten wurde der Wert jedoch erreicht.
RMSE erzielt einen Wert von 67195,66$ und ist damit auch eher hoch.

X_plot = X_plot.assign(pred_test2=lambda X_test: spline_pipe.predict(X_test))
# plot
sns.scatterplot(x=X_plot['median_income'], y=y_plot, label='actual')

plt.plot(X_plot['median_income'], X_plot['pred_test2'], label='predicted', color='orange', marker='.', linestyle='')
plt.legend();

Den Analysen entsprechend ist das beste Modell ist die Spline von Scikit Learn mit einem R-squared der Testdaten von 65,25%.
Gemäß dem Data-Science Lifecycle müsste dieses Modell jetzt deployed werden. Dieser Schritt wird allerdings nicht mehr im Rahmen dieses Projekts durchgeführt.

Applied Statistics Project

Regression

Contents

Regression¶

Python Setup¶

Datenaufbereitung¶

Daten laden¶

Daten inspizieren¶

Datentransformation¶

Entfernung der null values¶

Entfernung unwichtiger Features¶

Konvertierung von Datentypen¶

Hinzufügen neuer Features¶

Datensplit¶

Deskriptive Statistik¶

Kennzahlenanalyse numerischer Variablen¶

Analyse kategorialer Variablen¶

Exploratory Data Anaysis¶

erster Überblick¶

Korrelationsanalyse¶

Analyse der kategorialen Variable Ocean_Proximity¶

Analyse numerischer Variablen¶

Feature Auswahl¶

Modellierung¶

Statsmodel¶

Model 1 - OLS Regression¶

Regression Diagnostics¶

Model 2 - OLS Regression¶

Model 2 - Lasso Regression¶

Model 3 - Cubic Spline¶

Model 4 - Natural Spline¶

Scikit Learn¶

Data Preprocessing Pipeline¶

Datenvorbereitung¶

Model 5 - Simple Regression¶

Model 6 - Lasso¶

vorgegebenes Alpha¶

k-fold cross validation zur Ermittlung von alpha¶

Model 7 - Natural Splines¶

Evaluierung der besten 2 Modelle mit Test Datenset¶

Evaluation Statsmodels OLS lm2¶

Evaluation Scikit Learn Spline¶